#include <stdio.h>
#include <stdlib.h>
#include <string.h>

#include "galois.h"

#define NONE (10)
#define TABLE (11)
#define BITMATRIX (12)
#define LOGS (13)
#define SPLITW8 (14)

#define talloc(type, num) (type *) malloc(sizeof(type)*(num))

static int prim_poly[33] = 
{ 0, 
	/*  1 */     1, 
	/*  2 */    07,
	/*  3 */    013,
	/*  4 */    023,
	/*  5 */    045,
	/*  6 */    0103,
	/*  7 */    0211,
	/*  8 */    0435,
	/*  9 */    01021,
	/* 10 */    02011,
	/* 11 */    04005,
	/* 12 */    010123,
	/* 13 */    020033,
	/* 14 */    042103,
	/* 15 */    0100003,
	/* 16 */    0210013,
	/* 17 */    0400011,
	/* 18 */    01000201,
	/* 19 */    02000047,
	/* 20 */    04000011,
	/* 21 */    010000005,
	/* 22 */    020000003,
	/* 23 */    040000041,
	/* 24 */    0100000207,
	/* 25 */    0200000011,
	/* 26 */    0400000107,
	/* 27 */    01000000047,
	/* 28 */    02000000011,
	/* 29 */    04000000005,
	/* 30 */    010040000007,
	/* 31 */    020000000011, 
	/* 32 */    00020000007 };  /* Really 40020000007, but we're omitting the high order bit */

static int mult_type[33] = 
{ NONE, 
	/*  1 */   TABLE, 
	/*  2 */   TABLE,
	/*  3 */   TABLE,
	/*  4 */   TABLE,
	/*  5 */   TABLE,
	/*  6 */   TABLE,
	/*  7 */   TABLE,
	/*  8 */   TABLE,
	/*  9 */   TABLE,
	/* 10 */   LOGS,
	/* 11 */   LOGS,
	/* 12 */   LOGS,
	/* 13 */   LOGS,
	/* 14 */   LOGS,
	/* 15 */   LOGS,
	/* 16 */   LOGS,
	/* 17 */   LOGS,
	/* 18 */   LOGS,
	/* 19 */   LOGS,
	/* 20 */   LOGS,
	/* 21 */   LOGS,
	/* 22 */   LOGS,
	/* 23 */   BITMATRIX,
	/* 24 */   BITMATRIX,
	/* 25 */   BITMATRIX,
	/* 26 */   BITMATRIX,
	/* 27 */   BITMATRIX,
	/* 28 */   BITMATRIX,
	/* 29 */   BITMATRIX,
	/* 30 */   BITMATRIX,
	/* 31 */   BITMATRIX,
	/* 32 */   SPLITW8 };

static int nw[33] = { 0, (1 << 1), (1 << 2), (1 << 3), (1 << 4), 
	(1 << 5), (1 << 6), (1 << 7), (1 << 8), (1 << 9), (1 << 10),
	(1 << 11), (1 << 12), (1 << 13), (1 << 14), (1 << 15), (1 << 16),
	(1 << 17), (1 << 18), (1 << 19), (1 << 20), (1 << 21), (1 << 22),
	(1 << 23), (1 << 24), (1 << 25), (1 << 26), (1 << 27), (1 << 28),
	(1 << 29), (1 << 30), (1 << 31), -1 };

static int nwm1[33] = { 0, (1 << 1)-1, (1 << 2)-1, (1 << 3)-1, (1 << 4)-1, 
	(1 << 5)-1, (1 << 6)-1, (1 << 7)-1, (1 << 8)-1, (1 << 9)-1, (1 << 10)-1,
	(1 << 11)-1, (1 << 12)-1, (1 << 13)-1, (1 << 14)-1, (1 << 15)-1, (1 << 16)-1,
	(1 << 17)-1, (1 << 18)-1, (1 << 19)-1, (1 << 20)-1, (1 << 21)-1, (1 << 22)-1,
	(1 << 23)-1, (1 << 24)-1, (1 << 25)-1, (1 << 26)-1, (1 << 27)-1, (1 << 28)-1,
	(1 << 29)-1, (1 << 30)-1, 0x7fffffff, 0xffffffff };
	
static unsigned char *galois_w04_shifted_mult_table = { NULL };

static unsigned char *galois_w08_log_tables[9] = { NULL, NULL, NULL, NULL, NULL, NULL, NULL,
	NULL, NULL };

static uint16_t *galois_w16_log_tables[17] = { NULL, NULL, NULL, NULL, NULL, NULL, NULL,
	NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL };

static uint32_t *galois_w32_log_tables[33] = { NULL, NULL, NULL, NULL, NULL, NULL, NULL,
	NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL,
	NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL };

static unsigned char *galois_w08_ilog_tables[9] = { NULL, NULL, NULL, NULL, NULL, NULL, NULL,
	NULL, NULL };

static uint16_t *galois_w16_ilog_tables[17] = { NULL, NULL, NULL, NULL, NULL, NULL, NULL,
	NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL };

static uint32_t *galois_w32_ilog_tables[33] = { NULL, NULL, NULL, NULL, NULL, NULL, NULL,
	NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL,
	NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL };

static unsigned char *galois_w08_mult_tables[9] = { NULL, NULL, NULL, NULL, NULL, NULL, NULL,
	NULL, NULL };

static uint16_t *galois_w16_mult_tables[17] = { NULL, NULL, NULL, NULL, NULL, NULL, NULL,
	NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL };

static uint32_t *galois_w32_mult_tables[33] = { NULL, NULL, NULL, NULL, NULL, NULL, NULL,
	NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL,
	NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL };

static unsigned char *galois_w08_div_tables[9] = { NULL, NULL, NULL, NULL, NULL, NULL, NULL,
	NULL, NULL };

static uint16_t *galois_w16_div_tables[17] = { NULL, NULL, NULL, NULL, NULL, NULL, NULL,
	NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL };

static uint32_t *galois_w32_div_tables[33] = { NULL, NULL, NULL, NULL, NULL, NULL, NULL,
	NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL,
	NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL };

/* Special case for w = 32 */

static uint32_t *galois_w32_split8_tables[7] = { NULL, NULL, NULL, NULL, NULL, NULL, NULL };

static uint64_t *galois_w64_split8_tables[15] = { NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL };

static uint64_t *galois_w128_split8_tables[31] = { NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, 
												   NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL };

uint64_t rshift (uint64_t x, int shift) {int i; for (i = 0; i < shift; i++) x >>= 1; return x;}
uint64_t lshift (uint64_t x, int shift) {int i; for (i = 0; i < shift; i++) x <<= 1; return x;}

uint32_t galois_single_multiply(uint32_t x, uint32_t y, int w) {
	int sum_j;

	if (x == 0 || y == 0) return 0;

	if (mult_type[w] == TABLE) {
		if (w <= 8) {
			if (galois_w08_mult_tables[w] == NULL) {
				if (galois_w08_create_mult_tables(w) < 0) {
					fprintf(stderr, "ERROR -- cannot make multiplication tables for w=%d\n", w);
					exit(1);
				}
			}
			return (uint32_t) galois_w08_mult_tables[w][(x<<w)|y];

		} else {
			if (galois_w16_mult_tables[w] == NULL) {
				if (galois_w16_create_mult_tables(w) < 0) {
					fprintf(stderr, "ERROR -- cannot make multiplication tables for w=%d\n", w);
					exit(1);
				}
			}
			return (uint32_t) galois_w16_mult_tables[w][(x<<w)|y];
		}

	} else if (mult_type[w] == LOGS) {
		if (w <= 8) {
			if (galois_w08_log_tables[w] == NULL) {
				if (galois_w08_create_log_tables(w) < 0) {
					fprintf(stderr, "ERROR -- cannot make log tables for w=%d\n", w);
					exit(1);
				}
			}
			sum_j = galois_w08_log_tables[w][x] + galois_w08_log_tables[w][y];
			return (uint32_t) galois_w08_ilog_tables[w][sum_j];

		} else if (w <= 16) {
			if (galois_w16_log_tables[w] == NULL) {
				if (galois_w16_create_log_tables(w) < 0) {
					fprintf(stderr, "ERROR -- cannot make log tables for w=%d\n", w);
					exit(1);
				}
			}
			sum_j = galois_w16_log_tables[w][x] + galois_w16_log_tables[w][y];
			return (uint32_t) galois_w16_ilog_tables[w][sum_j];

		} else if (w <= 32) {
			if (galois_w32_log_tables[w] == NULL) {
				if (galois_w32_create_log_tables(w) < 0) {
					fprintf(stderr, "ERROR -- cannot make log tables for w=%d\n", w);
					exit(1);
				}
			}
			sum_j = galois_w32_log_tables[w][x] + galois_w32_log_tables[w][y];
			return (uint32_t) galois_w32_ilog_tables[w][sum_j];

		}
	} else if (mult_type[w] == SPLITW8) {
		if (galois_w32_split8_tables[0] == NULL) {
			if (galois_w32_create_split8_tables() < 0) {
				fprintf(stderr, "ERROR -- cannot make split8_tables for w=%d\n", w);
				exit(1);
			}
		}
		return galois_w32_split8_table_multiply(x, y);
	} else if (mult_type[w] == BITMATRIX) {
		return galois_bitmatrix_multiply(x, y, w);
	}
	fprintf(stderr, "Galois_single_multiply - no implementation for w=%d\n", w);
	exit(1);
}

uint32_t galois_single_divide(uint32_t a, uint32_t b, int w) {
	int sum_j;
	uint32_t inverse;

	if (mult_type[w] == TABLE) {
		if (w <= 8) {
			if (galois_w08_div_tables[w] == NULL) {
				if (galois_w08_create_mult_tables(w) < 0) {
					fprintf(stderr, "ERROR -- cannot make multiplication tables for w=%d\n", w);
					exit(1);
				}
			}
			return (uint32_t) galois_w08_div_tables[w][(a<<w)|b];

		} else if (w <= 16) {
			if (galois_w16_div_tables[w] == NULL) {
				if (galois_w16_create_mult_tables(w) < 0) {
					fprintf(stderr, "ERROR -- cannot make multiplication tables for w=%d\n", w);
					exit(1);
				}
			}
			return (uint32_t) galois_w16_div_tables[w][(a<<w)|b];

		}
	} else if (mult_type[w] == LOGS) {
		if (b == 0) return -1;
		if (a == 0) return 0;

		if (w <= 8) {
			if (galois_w08_log_tables[w] == NULL) {
				if (galois_w08_create_log_tables(w) < 0) {
					fprintf(stderr, "ERROR -- cannot make log tables for w=%d\n", w);
					exit(1);
				}
			}
			sum_j = galois_w08_log_tables[w][a] - galois_w08_log_tables[w][b];
			return (uint32_t) galois_w08_ilog_tables[w][sum_j];

		} else if (w <= 16) {
			if (galois_w16_log_tables[w] == NULL) {
				if (galois_w16_create_log_tables(w) < 0) {
					fprintf(stderr, "ERROR -- cannot make log tables for w=%d\n", w);
					exit(1);
				}
			}
			sum_j = galois_w16_log_tables[w][a] - galois_w16_log_tables[w][b];
			return (uint32_t) galois_w16_ilog_tables[w][sum_j];

		} else if (w <= 32) {
			if (galois_w32_log_tables[w] == NULL) {
				if (galois_w32_create_log_tables(w) < 0) {
					fprintf(stderr, "ERROR -- cannot make log tables for w=%d\n", w);
					exit(1);
				}
			}
			sum_j = galois_w32_log_tables[w][a] - galois_w32_log_tables[w][b];
			return (uint32_t) galois_w32_ilog_tables[w][sum_j];

		}
	} else {
		if (b == 0) return -1;
		if (a == 0) return 0;
		inverse = galois_inverse(b, w);
		return galois_single_multiply(a, inverse, w);
	}
	fprintf(stderr, "Galois_single_divide - no implementation for w=%d\n", w);
	exit(1);
}

uint32_t galois_inverse(uint32_t y, int w) {

	if (y == 0) return -1;
	if (mult_type[w] == BITMATRIX || mult_type[w] == SPLITW8) return galois_bitmatrix_inverse(y, w);
	return galois_single_divide(1, y, w);
}

int galois_create_log_tables(int w) {
	if (w > 30) return -1;
	if (w <= 8) {
		if (galois_w08_log_tables[w] == NULL) {
			if (galois_w08_create_log_tables(w) < 0) {
				fprintf(stderr, "Galois_w08_create_log_tables Error w: %d\n", w);
				return -1;
			}
		}
		return 0;
	} else if (w <= 16) {
		if (galois_w16_log_tables[w] == NULL) {
			if (galois_w16_create_log_tables(w) < 0) {
				fprintf(stderr, "Galois_w16_create_log_tables Error w: %d\n", w);
				return -1;
			}
		}
		return 0;
	} else if (w <= 32) {
		if (galois_w32_log_tables[w] == NULL) {
			if (galois_w32_create_log_tables(w) < 0) {
				fprintf(stderr, "Galois_w32_create_log_tables Error w: %d\n", w);
				return -1;
			}
		}
		return 0;
	} else {
		fprintf(stderr, "No logtable implementation for w=%d\n", w);
		return -1;
	}
}

uint32_t galois_logtable_multiply(uint32_t x, uint32_t y, int w) {
	int sum_j;

	if (x == 0 || y == 0) return 0;

	if (w <= 8) {
		sum_j = galois_w08_log_tables[w][x] + galois_w08_log_tables[w][y];
		return (uint32_t) galois_w08_ilog_tables[w][sum_j];
	} else if (w <= 16) {
		sum_j = galois_w16_log_tables[w][x] + galois_w16_log_tables[w][y];
		return (uint32_t) galois_w16_ilog_tables[w][sum_j];
	} else if (w <= 32) {
		sum_j = galois_w32_log_tables[w][x] + galois_w32_log_tables[w][y];
		return (uint32_t) galois_w32_ilog_tables[w][sum_j];
	}
}

uint32_t galois_logtable_divide(uint32_t x, uint32_t y, int w) {
	int sum_j;

	if (y == 0) return -1;
	if (x == 0) return 0;
	if (w <= 8) {
		sum_j = galois_w08_log_tables[w][x] - galois_w08_log_tables[w][y];
		return (uint32_t) galois_w08_ilog_tables[w][sum_j];

	} else if (w <= 16) {
		sum_j = galois_w16_log_tables[w][x] - galois_w16_log_tables[w][y];
		return (uint32_t) galois_w16_ilog_tables[w][sum_j];

	} else if (w <= 32) {
		sum_j = galois_w32_log_tables[w][x] - galois_w32_log_tables[w][y];
		return (uint32_t) galois_w32_ilog_tables[w][sum_j];
	}
}

int galois_create_mult_tables(int w) {
	if (w >= 14) return -1;
	if (w <= 8) {
		if (galois_w08_mult_tables[w] == NULL) {
			if (galois_w08_create_mult_tables(w) < 0) {
				fprintf(stderr, "Galois_w08_create_mult_tables Error w: %d\n", w);
				return -1;
			}
		}
		return 0;
	} else if (w <= 16) {
		if (galois_w16_mult_tables[w] == NULL) {
			if (galois_w16_create_mult_tables(w) < 0) {
				fprintf(stderr, "Galois_w16_create_mult_tables Error w: %d\n", w);
				return -1;
			}
		}
		return 0;
	} else {
		fprintf(stderr, "No multtable implementation for w=%d\n", w);
		return -1;
	}
}

uint32_t galois_multtable_multiply(uint32_t x, uint32_t y, int w) {
	if (w <= 8) {
		return (uint32_t) galois_w08_mult_tables[w][(x<<w)|y];
	} else if (w <= 16) {
		return (uint32_t) galois_w16_mult_tables[w][(x<<w)|y];
	} else {
		fprintf(stderr, "Galois_multtable_multiply - no implementation for w=%d\n", w);
		exit(1);
	}
}

uint32_t galois_multtable_divide(uint32_t x, uint32_t y, int w) {
	if (w <= 8) {
		if (y == 0) return -1;
		return (uint32_t) galois_w08_div_tables[w][(x<<w)|y];
	} else if (w <= 16) {
		return (uint32_t) galois_w16_div_tables[w][(x<<w)|y];
	} else {
		fprintf(stderr, "Galois_multtable_divide - no implementation for w=%d\n", w);
		exit(1);
	}
}

uint32_t galois_bitmatrix_multiply(uint32_t x, uint32_t y, int w) {
	uint32_t prod;
	int i, ind;
	int k;
	uint32_t scratch[32];

	prod = 0;
	for (i = 0; i < w; i++) {
		scratch[i] = y;
		if (y & (1 << (w-1))) {
			y = y << 1;
			y = (y ^ prim_poly[w]) & nwm1[w];
		} else {
			y = y << 1;
		}
	}
	for (i = 0; i < w; i++) {
		ind = (1 << i);
		if (ind & x) {
			prod = prod ^ scratch[i];
		}
	}
	return prod;
}

uint32_t galois_bitmatrix_divide(uint32_t a, uint32_t b, int w) {
	uint32_t inverse;

	if (b == 0) return -1;
	if (a == 0) return 0;
	inverse = galois_bitmatrix_inverse(b, w);
	return galois_bitmatrix_multiply(a, inverse, w);
}

/* This will destroy mat, by the way */
void galois_invert_binary_matrix(int *mat, int *inv, int rows) {
	int cols, i, j, k;
	int tmp;

	cols = rows;

	for (i = 0; i < rows; i++) inv[i] = (1 << i);

	/* First -- convert into upper triangular */

	for (i = 0; i < cols; i++) {

		/* Swap rows if we ave a zero i,i element.  If we can't swap, then the 
			 matrix was not invertible */

		if ((mat[i] & (1 << i)) == 0) { 
			for (j = i+1; j < rows && (mat[j] & (1 << i)) == 0; j++) ;
			if (j == rows) {
				fprintf(stderr, "galois_invert_matrix: Matrix not invertible!!\n");
				exit(1);
			}
			tmp = mat[i]; mat[i] = mat[j]; mat[j] = tmp;
			tmp = inv[i]; inv[i] = inv[j]; inv[j] = tmp;
		}

		/* Now for each j>i, add A_ji*Ai to Aj */
		for (j = i+1; j != rows; j++) {
			if ((mat[j] & (1 << i)) != 0) {
				mat[j] ^= mat[i]; 
				inv[j] ^= inv[i];
			}
		}
	}

	/* Now the matrix is upper triangular.  Start at the top and multiply down */

	for (i = rows-1; i >= 0; i--) {
		for (j = 0; j < i; j++) {
			if (mat[j] & (1 << i)) {
				/*        mat[j] ^= mat[i]; */
				inv[j] ^= inv[i];
			}
		}
	} 
}

uint32_t galois_bitmatrix_inverse(uint32_t y, int w) {
	int mat[1024], mat2[32];
	int inv[1024], inv2[32];
	int ind, i, j, k, prod;

	for (i = 0; i < w; i++) {
		mat2[i] = y;

		if (y & nw[w-1]) {
			y = y << 1;
			y = (y ^ prim_poly[w]) & nwm1[w];
		} else {
			y = y << 1;
		}
	}

	galois_invert_binary_matrix(mat2, inv2, w);

	return (uint32_t) inv2[0]; 
}

uint32_t galois_shift_multiply(uint32_t x, uint32_t y, int w) {
	uint64_t product, one, n2, pp;
	int i;

	n2 = y;
	product = 0;
	one = 1;
	pp = prim_poly[w];
	if (w == 32) pp |= (one << 32);
	for (i = 0; i < w; i++) {
		if ((x & (1 << i))) product ^= (n2 << i);
	}
	for (i = 2*w - 1; i >= w; i--) {
		if ((product & (one << i))) product ^= (pp << (i-w));
	}
	i = product & 0xffffffff;
	return i;
}

uint32_t galois_shift_divide (uint32_t x, uint32_t y, int w) {
	uint32_t inverse;

	if (y == 0) return -1;
	if (x == 0) return 0;
	inverse = galois_euclid_inverse(y, w);
	return galois_shift_multiply(x, inverse, w);

}

uint32_t galois_euclid_inverse(uint32_t x, int w) {
	int remainder[2], auxiliary[2];

	int w32;
	int high_bit;
	int high_bit_r;
	int r1, q1, shift_count, a1;

	if (x == 1) return 1;
	
	remainder[0] = prim_poly[w];
	remainder[1] = x;

	auxiliary[0] = 0;
	auxiliary[1] = 1;
	
	w32 = (w == 32) ? 1 : 0;

	high_bit_r = w;
	while (remainder[1] !=1) {
		q1 = 0;
		high_bit = high_bit_r;
		if (w32 == 1) {
			high_bit_r--;
		}

		while((remainder[1] & (1 << high_bit_r)) == 0) high_bit_r--;

		/* This is the polynomial division part */
		while(high_bit >= high_bit_r) {
			r1 = remainder[1];

			shift_count = high_bit - high_bit_r;
			r1 <<= shift_count;
			

			/* set the quotient */
			q1 ^= (1 << shift_count);

			/* next remainder in the long division process */
			remainder[0] ^= r1;
			if (w32 == 1) { high_bit--; w32 = 0; }
			while((remainder[0] & (1 << high_bit)) == 0) high_bit--;
		}

		r1 = remainder[0];
		remainder[0] = remainder[1];
		remainder[1] = r1;
		a1 = galois_shift_multiply(q1, auxiliary[1], w);
		a1 ^= auxiliary[0];
		auxiliary[0] = auxiliary[1];
		auxiliary[1] = a1;
		
	}
	
	return a1;
}

uint32_t galois_ilog(uint32_t value, int w) {
	if (w <= 8) {
		if (galois_w08_ilog_tables[w] == NULL) {
			if (galois_w08_create_log_tables(w) < 0) {
				fprintf(stderr, "Error: galois_ilog - w is too big.  Sorry\n");
				exit(1);
			}
		}
		return (uint32_t) galois_w08_ilog_tables[w][value];

	} else if (w <= 16) {
		if (galois_w16_ilog_tables[w] == NULL) {
			if (galois_w16_create_log_tables(w) < 0) {
				fprintf(stderr, "Error: galois_ilog - w is too big.  Sorry\n");
				exit(1);
			}
		}
		return (uint32_t) galois_w16_ilog_tables[w][value];

	} else if (w <= 32) {
		if (galois_w32_ilog_tables[w] == NULL) {
			if (galois_w32_create_log_tables(w) < 0) {
				fprintf(stderr, "Error: galois_ilog - w is too big.  Sorry\n");
				exit(1);
			}
		}
		return (uint32_t) galois_w32_ilog_tables[w][value];

	}
	fprintf(stderr, "Galois_ilog - no implementation for w=%d\n", w);
	exit(1);
}

uint32_t galois_log(uint32_t value, int w) {
	if (w <= 8) {
		if (galois_w08_log_tables[w] == NULL) {
			if (galois_w08_create_log_tables(w) < 0) {
				fprintf(stderr, "Error: galois_log - w is too big.  Sorry\n");
				exit(1);
			}
		}
		return (uint32_t) galois_w08_log_tables[w][value];

	} else if (w <= 16) {
		if (galois_w16_log_tables[w] == NULL) {
			if (galois_w16_create_log_tables(w) < 0) {
				fprintf(stderr, "Error: galois_log - w is too big.  Sorry\n");
				exit(1);
			}
		}
		return (uint32_t) galois_w16_log_tables[w][value];

	} else if (w <= 32) {
		if (galois_w32_log_tables[w] == NULL) {
			if (galois_w32_create_log_tables(w) < 0) {
				fprintf(stderr, "Error: galois_log - w is too big.  Sorry\n");
				exit(1);
			}
		}
		return (uint32_t) galois_w32_log_tables[w][value];

	}
	fprintf(stderr, "Galois_log - no implementation for w=%d\n", w);
	exit(1);
}

/* Single operations for minimum width functions */
/* w <= 8 */

int galois_w08_create_log_tables (int w) {
	int j, b;

	if (w > 8) return -1;
	if (galois_w08_log_tables[w] != NULL) return 0;
	galois_w08_log_tables[w] = (unsigned char *) malloc(sizeof(unsigned char)*nw[w]);
	if (galois_w08_log_tables[w] == NULL) return -1; 

	galois_w08_ilog_tables[w] = (unsigned char *) malloc(sizeof(unsigned char)*nw[w]*3);
	if (galois_w08_ilog_tables[w] == NULL) { 
		free(galois_w08_log_tables[w]);
		galois_w08_log_tables[w] = NULL;
		return -1;
	}

	for (j = 0; j < nw[w]; j++) {
		galois_w08_log_tables[w][j] = (unsigned char) nwm1[w];
		galois_w08_ilog_tables[w][j] = (unsigned char) 0;
	} 

	b = 1;
	for (j = 0; j < nwm1[w]; j++) {
		if (galois_w08_log_tables[w][b] != (unsigned char) nwm1[w]) {
			fprintf(stderr, "Galois_create_w08_log_tables Error: j=%d, b=%d, B->J[b]=%d, J->B[j]=%d (0%o)\n",
					j, b, (int) galois_w08_log_tables[w][b], (int) galois_w08_ilog_tables[w][j], (b << 1) ^ prim_poly[w]);
			exit(1);
		}
		galois_w08_log_tables[w][b] = (unsigned char) j;
		galois_w08_ilog_tables[w][j] = (unsigned char) b;
		b = b << 1;
		if (b & nw[w]) b = (b ^ prim_poly[w]) & nwm1[w];
	}
	for (j = 0; j < nwm1[w]; j++) {
		galois_w08_ilog_tables[w][j+nwm1[w]] = galois_w08_ilog_tables[w][j];
		galois_w08_ilog_tables[w][j+nwm1[w]*2] = galois_w08_ilog_tables[w][j];
	} 
	galois_w08_ilog_tables[w] += nwm1[w];
	return 0;

}

unsigned char galois_w08_logtable_multiply(unsigned char x, unsigned char y, int w) {
	int sum_j;

	if (x == 0 || y == 0) return 0;

	sum_j = galois_w08_log_tables[w][x] + galois_w08_log_tables[w][y];
	/* if (sum_j >= nwm1[w]) sum_j -= nwm1[w];    Don't need to do this, 
		 because we replicate the ilog table twice.  */
	return galois_w08_ilog_tables[w][sum_j];
}

unsigned char galois_w08_logtable_divide(unsigned char x, unsigned char y, int w) {
	int sum_j;
	unsigned char z;

	if (y == 0) return -1;
	if (x == 0) return 0; 
	sum_j = galois_w08_log_tables[w][x] - galois_w08_log_tables[w][y];
	/* if (sum_j < 0) sum_j += nwm1[w];   Don't need to do this, because we replicate the ilog table twice.   */
	z = galois_w08_ilog_tables[w][sum_j];
	return z;
}

int galois_w08_create_mult_tables(int w) {
	int j, x, y, logx;

	if (w > 8) return -1;

	if (galois_w08_mult_tables[w] != NULL) return 0;
	galois_w08_mult_tables[w] = (unsigned char *) malloc(sizeof(unsigned char) * nw[w] * nw[w]);
	if (galois_w08_mult_tables[w] == NULL) return -1;

	galois_w08_div_tables[w] = (unsigned char *) malloc(sizeof(unsigned char) * nw[w] * nw[w]);
	if (galois_w08_div_tables[w] == NULL) {
		free(galois_w08_mult_tables[w]);
		galois_w08_mult_tables[w] = NULL;
		return -1;
	}
	if (galois_w08_log_tables[w] == NULL) {
		if (galois_w08_create_log_tables(w) < 0) {
			free(galois_w08_mult_tables[w]);
			free(galois_w08_div_tables[w]);
			galois_w08_mult_tables[w] = NULL;
			galois_w08_div_tables[w] = NULL;
			return -1;
		}
	}

	/* Set mult/div tables for x = 0 */
	j = 0;
	galois_w08_mult_tables[w][j] = (unsigned char) 0;   /* y = 0 */
	galois_w08_div_tables[w][j] = (unsigned char) -1;
	j++;
	for (y = 1; y < nw[w]; y++) {   /* y > 0 */
		galois_w08_mult_tables[w][j] = (unsigned char) 0;
		galois_w08_div_tables[w][j] = (unsigned char) 0;
		j++;
	}

	for (x = 1; x < nw[w]; x++) {  /* x > 0 */
		galois_w08_mult_tables[w][j] = (unsigned char) 0; /* y = 0 */
		galois_w08_div_tables[w][j] = (unsigned char) -1;
		j++;
		logx = (int) galois_w08_log_tables[w][x];
		for (y = 1; y < nw[w]; y++) {  /* y > 0 */
			galois_w08_mult_tables[w][j] = (unsigned char) galois_w08_ilog_tables[w][logx+(int) galois_w08_log_tables[w][y]]; 
			galois_w08_div_tables[w][j]  = (unsigned char) galois_w08_ilog_tables[w][logx-(int) galois_w08_log_tables[w][y]]; 
			j++;
		}
	}
	return 0;
}

unsigned char galois_w08_multtable_multiply (unsigned char x, unsigned char y, int w) {
	return galois_w08_mult_tables[w][(x<<w)|y];
}

int galois_w08_multtable_divide (unsigned char x, unsigned char y, int w) {
	if (y == 0) return -1;
	return (int) galois_w08_div_tables[w][(x<<w)|y];
}

unsigned char galois_w08_composite_k2_multiply(unsigned char x, unsigned char y) {
	if (x == 0 || y == 0) return 0;
	int l;
	unsigned char s;
	unsigned char product, mask0;
	unsigned char p2, p3;
	
	unsigned char *table;

	if (galois_w08_mult_tables[4] == NULL) {
		if (galois_w08_create_mult_tables(4) < 0) {
			fprintf(stderr, "galois_w08_composite_k2_multiply: Error creating tables for w = 4\n");
			exit(1);
		}
	}

	l = 4;
	s = 4;
	mask0 = (1 << l) - 1;
	table = galois_w08_mult_tables[4];

	p2 = table[((x>>l)<<l)|(y>>l)];
	p3 = p2;
	p2 = table[(s<<l)|p2];

	product = table[((x>>l)<<l)|(y & mask0)];
	product ^= table[((x & mask0) << l) | (y>>l)];
	product ^= p2;
	product <<= l;
	product ^= p3;
	product ^= table[((x & mask0) << l) | (y & mask0)];

	return product;
}

unsigned char galois_w08_composite_k2_divide(unsigned char x, unsigned char y) {
	unsigned char inverse;
	if (y == 0) return -1;
	if (x == 0) return 0;

	inverse = galois_w08_composite_k2_inverse (y);
	return galois_w08_composite_k2_multiply (x, inverse);

}

unsigned char galois_w08_composite_k2_inverse(unsigned char x) {
	if (x == 0) return -1;
	int l;
	unsigned char s;
	unsigned char c;
	unsigned char mask0, inv;
	unsigned char a1i, a0i;

	if (galois_w08_div_tables[4] == NULL) {
		if (galois_w08_create_mult_tables(4) < 0) {
			fprintf(stderr, "galois_w08_composite_k2_inverse: Error creating mult tables for w=4\n");
			exit(1);
		}
	}

	l = 4;
	s = 4;
	mask0 = (1 << l) - 1;

	if ((x >> l) & mask0 == 0) {
		inv = galois_w08_div_tables[4][(1<<l) | (x & mask0)];
	} else if ((x & mask0) == 0) {
		a1i = galois_w08_div_tables[4][(1<<l) | ((x >> l) & mask0)];
		inv = a1i;
		inv <<= l;
		inv ^= galois_w08_mult_tables[4][(a1i<<l) | s];
	} else {
		a1i = galois_w08_div_tables[4][(1 << l) | ((x >> l) & mask0)];
		a0i = galois_w08_div_tables[4][(1 << l) | (x & mask0)];
		
		c = galois_w08_mult_tables[4][(a0i << l) | ((x >> l) &  mask0)];
		c ^= galois_w08_mult_tables[4][((x & mask0) << l) | a1i];
		c ^= s & mask0;
		c = galois_w08_div_tables[4][(1 << l) | (c & mask0)];
		c = galois_w08_mult_tables[4][((c & mask0) << l) | ((x >> l) & mask0)];
		c = galois_w08_mult_tables[4][(a0i << l) | (c & mask0)];

		inv = galois_w08_mult_tables[4][(a1i << l) | (c & mask0)];
		inv <<= l;
		inv ^= galois_w08_mult_tables[4][(a0i << l) | (c ^ 1)];

	}
	return inv;

}

/* w <= 16 */

int galois_w16_create_log_tables (int w) {
	int j, b;

	if (w > 16) return -1;
	if (galois_w16_log_tables[w] != NULL) return 0;
	galois_w16_log_tables[w] = (uint16_t *) malloc(sizeof(uint16_t)*nw[w]);
	if (galois_w16_log_tables[w] == NULL) return -1; 

	galois_w16_ilog_tables[w] = (uint16_t *) malloc(sizeof(uint16_t)*nw[w]*3);
	if (galois_w16_ilog_tables[w] == NULL) { 
		free(galois_w16_log_tables[w]);
		galois_w16_log_tables[w] = NULL;
		return -1;
	}

	for (j = 0; j < nw[w]; j++) {
		galois_w16_log_tables[w][j] = (uint16_t) nwm1[w];
		galois_w16_ilog_tables[w][j] = (uint16_t) 0;
	} 

	b = 1;
	for (j = 0; j < nwm1[w]; j++) {
		if (galois_w16_log_tables[w][b] != (uint16_t) nwm1[w]) {
			fprintf(stderr, "Galois_create_w16_log_tables Error: j=%d, b=%d, B->J[b]=%d, J->B[j]=%d (0%o)\n",
					j, b, (int) galois_w16_log_tables[w][b], (int) galois_w16_ilog_tables[w][j], (b << 1) ^ prim_poly[w]);
			exit(1);
		}
		galois_w16_log_tables[w][b] = (uint16_t) j;
		galois_w16_ilog_tables[w][j] = (uint16_t) b;
		b = b << 1;
		if (b & nw[w]) b = (b ^ prim_poly[w]) & nwm1[w];
	}
	for (j = 0; j < nwm1[w]; j++) {
		galois_w16_ilog_tables[w][j+nwm1[w]] = galois_w16_ilog_tables[w][j];
		galois_w16_ilog_tables[w][j+nwm1[w]*2] = galois_w16_ilog_tables[w][j];
	} 
	galois_w16_ilog_tables[w] += nwm1[w];
	return 0;

}

uint16_t galois_w16_logtable_multiply(uint16_t x, uint16_t y, int w) {
	int sum_j;

	if (x == 0 || y == 0) return 0;

	sum_j = galois_w16_log_tables[w][x] + galois_w16_log_tables[w][y];
	return galois_w16_ilog_tables[w][sum_j];
}

uint16_t galois_w16_logtable_divide(uint16_t x, uint16_t y, int w) {
	int sum_j;
	uint16_t z;

	if (y == 0) return -1;
	if (x == 0) return 0; 
	sum_j = galois_w16_log_tables[w][x] - galois_w16_log_tables[w][y];
	z = galois_w16_ilog_tables[w][sum_j];
	return z;
}

int galois_w16_create_mult_tables(int w) {
	int j, x, y, logx;

	if (w > 16) return -1;

	if (galois_w16_mult_tables[w] != NULL) return 0;
	galois_w16_mult_tables[w] = (uint16_t *) malloc(sizeof(uint16_t) * nw[w] * nw[w]);
	if (galois_w16_mult_tables[w] == NULL) return -1;

	galois_w16_div_tables[w] = (uint16_t *) malloc(sizeof(uint16_t) * nw[w] * nw[w]);
	if (galois_w16_div_tables[w] == NULL) {
		free(galois_w16_mult_tables[w]);
		galois_w16_mult_tables[w] = NULL;
		return -1;
	}
	if (galois_w16_log_tables[w] == NULL) {
		if (galois_w16_create_log_tables(w) < 0) {
			free(galois_w16_mult_tables[w]);
			free(galois_w16_div_tables[w]);
			galois_w16_mult_tables[w] = NULL;
			galois_w16_div_tables[w] = NULL;
			return -1;
		}
	}

	/* Set mult/div tables for x = 0 */
	j = 0;
	galois_w16_mult_tables[w][j] = (uint16_t) 0;   /* y = 0 */
	galois_w16_div_tables[w][j] = (uint16_t) -1;
	j++;
	for (y = 1; y < nw[w]; y++) {   /* y > 0 */
		galois_w16_mult_tables[w][j] = (uint16_t) 0;
		galois_w16_div_tables[w][j] = (uint16_t) 0;
		j++;
	}

	for (x = 1; x < nw[w]; x++) {  /* x > 0 */
		galois_w16_mult_tables[w][j] = (uint16_t) 0; /* y = 0 */
		galois_w16_div_tables[w][j] = (uint16_t) -1;
		j++;
		logx = (int) galois_w16_log_tables[w][x];
		for (y = 1; y < nw[w]; y++) {  /* y > 0 */
			galois_w16_mult_tables[w][j] = (uint16_t) galois_w16_ilog_tables[w][logx+(int) galois_w16_log_tables[w][y]]; 
			galois_w16_div_tables[w][j]  = (uint16_t) galois_w16_ilog_tables[w][logx-(int) galois_w16_log_tables[w][y]]; 
			j++;
		}
	}
	return 0;
}

uint16_t galois_w16_multtable_multiply (uint16_t x, uint16_t y, int w) {
	return galois_w16_mult_tables[w][(x<<w)|y];
}

uint16_t galois_w16_multtable_divide (uint16_t x, uint16_t y, int w) {
	return galois_w16_mult_tables[w][(w<<w)|y];
}

uint16_t galois_w16_composite_k2_multiply(uint16_t x, uint16_t y) {
	if (x == 0 || y == 0) return 0;
	int l;
	uint16_t s;
	uint16_t product, mask0;
	uint16_t p2, p3;
	
	unsigned char *table;

	if (galois_w08_mult_tables[8] == NULL) {
		if (galois_w08_create_mult_tables(8) < 0) {
			fprintf(stderr, "galois_w16_composite_k2_multiply: Error creating tables for w = 8\n");
			exit(1);
		}
	}

	l = 8;
	s = 63;
	mask0 = (1 << l) - 1;
	table = galois_w08_mult_tables[8];

	p2 = table[((x>>l)<<l)|(y>>l)];
	p3 = p2;
	p2 = table[(s<<l)|p2];

	product = table[((x>>l)<<l)|(y & mask0)];
	product ^= table[((x & mask0) << l) | (y>>l)];
	product ^= p2;
	product <<= l;
	product ^= p3;
	product ^= table[((x & mask0) << l) | (y & mask0)];

	return product;

}

uint16_t galois_w16_composite_k2_divide(uint16_t x, uint16_t y) {
	uint16_t inverse;
	if (y == 0) return -1;
	if (x == 0) return 0;

	inverse = galois_w16_composite_k2_inverse (y);
	return galois_w16_composite_k2_multiply (x, inverse);
	
}

uint16_t galois_w16_composite_k2_inverse(uint16_t x) {
	if (x == 0) return -1;
	int l;
	uint16_t s;
	uint16_t c;
	uint16_t mask0, inv;
	uint16_t a1i, a0i;

	if (galois_w08_div_tables[8] == NULL) {
		if (galois_w08_create_mult_tables(8) < 0) {
			fprintf(stderr, "galois_w08_composite_k2_inverse: Error creating mult tables for w=8\n");
			exit(1);
		}
	}

	l = 8;
	s = 63;
	mask0 = (1 << l) - 1;

	if ((x >> l) & mask0 == 0) {
		inv = galois_w08_div_tables[8][(1<<l) | (x & mask0)];
	} else if ((x & mask0) == 0) {
		a1i = galois_w08_div_tables[8][(1<<l) | ((x >> l) & mask0)];
		inv = a1i;
		inv <<= l;
		inv ^= galois_w08_mult_tables[8][(a1i<<l) | s];
	} else {
		a1i = galois_w08_div_tables[8][(1 << l) | ((x >> l) & mask0)];
		a0i = galois_w08_div_tables[8][(1 << l) | (x & mask0)];
		
		c = galois_w08_mult_tables[8][(a0i << l) | ((x >> l) &  mask0)];
		c ^= galois_w08_mult_tables[8][((x & mask0) << l) | a1i];
		c ^= s & mask0;
		c = galois_w08_div_tables[8][(1 << l) | (c & mask0)];
		c = galois_w08_mult_tables[8][((c & mask0) << l) | ((x >> l) & mask0)];
		c = galois_w08_mult_tables[8][(a0i << l) | (c & mask0)];

		inv = galois_w08_mult_tables[8][(a1i << l) | (c & mask0)];
		inv <<= l;
		inv ^= galois_w08_mult_tables[8][(a0i << l) | (c ^ 1)];

	}
	return inv;
	
}

/* w <= 32 */

int galois_w32_create_log_tables (int w) {
	int j, b;

	if (w > 32) return -1;
	if (galois_w32_log_tables[w] != NULL) return 0;
	galois_w32_log_tables[w] = (uint32_t *) malloc(sizeof(uint32_t)*nw[w]);
	if (galois_w32_log_tables[w] == NULL) return -1; 

	galois_w32_ilog_tables[w] = (uint32_t *) malloc(sizeof(uint32_t)*nw[w]*3);
	if (galois_w32_ilog_tables[w] == NULL) { 
		free(galois_w32_log_tables[w]);
		galois_w32_log_tables[w] = NULL;
		return -1;
	}

	for (j = 0; j < nw[w]; j++) {
		galois_w32_log_tables[w][j] = (uint32_t) nwm1[w];
		galois_w32_ilog_tables[w][j] = (uint32_t) 0;
	} 

	b = 1;
	for (j = 0; j < nwm1[w]; j++) {
		if (galois_w32_log_tables[w][b] != (uint32_t) nwm1[w]) {
			fprintf(stderr, "Galois_create_w32_log_tables Error: j=%d, b=%d, B->J[b]=%d, J->B[j]=%d (0%o)\n",
					j, b, (int) galois_w32_log_tables[w][b], (int) galois_w32_ilog_tables[w][j], (b << 1) ^ prim_poly[w]);
			exit(1);
		}
		galois_w32_log_tables[w][b] = (uint32_t) j;
		galois_w32_ilog_tables[w][j] = (uint32_t) b;
		b = b << 1;
		if (b & nw[w]) b = (b ^ prim_poly[w]) & nwm1[w];
	}
	for (j = 0; j < nwm1[w]; j++) {
		galois_w32_ilog_tables[w][j+nwm1[w]] = galois_w32_ilog_tables[w][j];
		galois_w32_ilog_tables[w][j+nwm1[w]*2] = galois_w32_ilog_tables[w][j];
	} 
	galois_w32_ilog_tables[w] += nwm1[w];
	return 0;

}

uint32_t galois_w32_logtable_multiply(uint32_t x, uint32_t y, int w) {
	int sum_j;

	if (x == 0 || y == 0) return 0;

	sum_j = galois_w32_log_tables[w][x] + galois_w32_log_tables[w][y];
	return galois_w32_ilog_tables[w][sum_j];
}

uint32_t galois_w32_logtable_divide(uint32_t x, uint32_t y, int w) {
	int sum_j;
	uint32_t z;

	if (y == 0) return -1;
	if (x == 0) return 0; 
	sum_j = galois_w32_log_tables[w][x] - galois_w32_log_tables[w][y];
	z = galois_w32_ilog_tables[w][sum_j];
	return z;
}

int galois_w32_create_split8_tables() {
	int p1, p2, i, j, index, ishift, jshift;
	uint32_t p1elt, p2elt, *table;

	if (galois_w32_split8_tables[0] != NULL) return 0;

	for (i = 0; i < 7; i++) {
		galois_w32_split8_tables[i] = (uint32_t *) malloc(sizeof(uint32_t) * (1 << 16));
		if (galois_w32_split8_tables[i] == NULL) {
			for (i--; i >= 0; i--) free(galois_w32_split8_tables[i]);
			return -1;
		}
	}

	for (i = 0; i < 4; i += 3) {
		ishift = i * 8;
		for (j = ((i == 0) ? 0 : 1) ; j < 4; j++) {
			jshift = j * 8;
			table = galois_w32_split8_tables[i+j];
			index = 0;
			for (p1 = 0; p1 < 256; p1++) {
				p1elt = (p1 << ishift);
				for (p2 = 0; p2 < 256; p2++) {
					p2elt = (p2 << jshift);
					table[index] = galois_bitmatrix_multiply(p1elt, p2elt, 32);
					index++;
				}
			}
		}
	}
	return 0;
}

uint32_t galois_w32_split8_table_multiply(uint32_t x, uint32_t y) {
	int i, j, a, b, i8, j8;
	uint32_t accumulator;

	accumulator = 0;

	i8 = 0;
	for (i = 0; i < 4; i++) {
		a = (((x >> i8) & 255) << 8);
		j8 = 0;
		for (j = 0; j < 4; j++) {
			b = ((y >> j8) & 255);
			accumulator ^= galois_w32_split8_tables[i+j][a|b];
			j8 += 8;
		}
		i8 += 8;
	}
	return accumulator;
}

uint32_t galois_w32_composite_k2_multiply(uint32_t x, uint32_t y) {
	if (x == 0 || y == 0) return 0;
	int l;
	uint32_t s;
	uint32_t product, mask0;
	uint32_t p2, p3;
	
	uint16_t *log_table, *ilog_table;

	if (galois_w16_log_tables[16] == NULL) {
		if (galois_w16_create_log_tables(16) < 0) {
			fprintf(stderr, "galois_w32_composite_k2_multiply: Error creating tables for w = 16\n");
			exit(1);
		}
	}

	l = 16;
	s = 40188;
	mask0 = (1 << l) - 1;
	log_table = galois_w16_log_tables[16];
	ilog_table = galois_w16_ilog_tables[16];

	p2 = galois_w16_logtable_multiply((x >> l), (y >> l), l);
	p3 = p2;
	p2 = galois_w16_logtable_multiply(s, p2, l);
	
	product = galois_w16_logtable_multiply((x >> l), (y & mask0), l);
	product ^= galois_w16_logtable_multiply((x & mask0), (y >> l), l);
	product ^= p2;
	product <<= l;
	product ^= p3;
	product ^= galois_w16_logtable_multiply((x & mask0), (y & mask0), l);
	
	/*
	p2 = ilog_table[ log_table[(x >> l)] + log_table[(y >> l)] ];
	p3 = p2;
	p2 = ilog_table[ log_table[s] + log_table[p2] ];

	product  = ilog_table[ log_table[(x >> l)]    + log_table[(y & mask0)] ];
	product ^= ilog_table[ log_table[(x & mask0)] + log_table[(y>>l)] ];
	product ^= p2;
	product <<= l;
	product ^= p3;
	product ^= ilog_table[ log_table[(x & mask0)] + log_table[(y & mask0)] ];
	*/

	return product;
	
}

uint32_t galois_w32_composite_k2_divide(uint32_t x, uint32_t y) {
	uint32_t inverse;
	if (y == 0) return -1;
	if (x == 0) return 0;

	inverse = galois_w32_composite_k2_inverse (y);
	return galois_w32_composite_k2_multiply (x, inverse);
	
}

uint32_t galois_w32_composite_k2_inverse(uint32_t x) {
	if (x == 0) return -1;
	int l;
	uint32_t s;
	uint32_t c;
	uint32_t mask0, inv;
	uint32_t a1i, a0i;
	int sum_log;

	if (galois_w16_log_tables[16] == NULL) {
		if (galois_w16_create_log_tables(16) < 0) {
			fprintf(stderr, "galois_w32_composite_k2_inverse: Error creating mult tables for w=16\n");
			exit(1);
		}
	}

	l = 16;
	s = 40188;
	mask0 = (1 << l) - 1;

	if ((x >> l) & mask0 == 0) {
		inv = galois_w16_logtable_divide(1, (x >> l), l);
		
	} else if ((x & mask0) == 0) {
		a1i = galois_w16_logtable_divide(1, (x >> l), l);
		inv = a1i;
		inv <<= l;
		inv ^= galois_w16_logtable_multiply(a1i, s, l);
		
	} else {
		a1i = galois_w16_logtable_divide(1, ((x>>l) & mask0), l);
		a0i = galois_w16_logtable_divide(1, (x &mask0), l);
		
		c = galois_w16_logtable_multiply(a0i, ((x >> l) & mask0), l);
		c^= galois_w16_logtable_multiply((x & mask0), a1i, l);
		c^= s & mask0;
		c = galois_w16_logtable_divide(1, c, l);
		c = galois_w16_logtable_multiply((c & mask0), ((x >> l) & mask0), l);
		c = galois_w16_logtable_multiply(a0i, (c & mask0), l);
		
		inv = galois_w16_logtable_multiply(a1i, (c & mask0), l);
		inv <<= l;
		inv^= galois_w16_logtable_multiply(a0i, c^1, l);

	}
	return inv;
	
}

/* w == 64 */

uint64_t galois_w64_bitmatrix_multiply (uint64_t x, uint64_t y) {
	int i, k;
	uint64_t prod, j;
	uint64_t scratch[64];
	uint64_t one;
	uint64_t pp;
	uint64_t wshift;

	one = 1;
	wshift = lshift(one, 63);
	pp = 0;
	pp = 0x1b;

	for (i = 0; i < 64; i++) {
		scratch[i] = y;

		if ((y & wshift)) {
			y = y << 1;
			y ^= pp;
		} else {
			y = y << 1;
		}
	}

	prod = 0;
	for (i = 0; i < 64; i++) {
		if ((x & (one << i))) {			
			prod ^= scratch[i];
		}
	}

	return prod;
}

uint64_t galois_w64_bitmatrix_divide (uint64_t x, uint64_t y) {
	uint64_t inverse;
	if (y == 0) return -1;
	if (x == 0) return 0;

	inverse = galois_w64_bitmatrix_inverse(y);
	return galois_w64_bitmatrix_multiply(x, inverse);
}

void galois_w64_invert_binary_matrix (uint64_t *mat, uint64_t *inv) {
	uint64_t cols, rows;
	int i, j, k;
	uint64_t tmp;
	uint64_t one;

	one = 1;
	cols = 64;
	rows = 64;

	for (i = 0; i < rows; i++) {
		inv[i] = 0;
		inv[i] = (one << i);
	}

	/* First -- convert into upper triangular */

	for (i = 0; i < cols; i++) {

		/* Swap rows if we ave a zero i,i element.  If we can't swap, then the 
			 matrix was not invertible */

		if ((mat[i] & (one << i)) == 0) {
			for (j = i+1; 
					j < rows && (mat[j] & (one << i)) == 0; 
					j++);
			if (j == rows) {
				fprintf(stderr, "galois_invert_matrix: Matrix not invertible!!\n");
				exit(1);
			}
			tmp = mat[i]; mat[i] = mat[j]; mat[j] = tmp;
			tmp = inv[i]; inv[i] = inv[j]; inv[j] = tmp;
		}

		/* Now for each j>i, add A_ji*Ai to Aj */
		for (j = i+1; j != rows; j++) {
			if ((mat[j] & (one << i)) != 0) {
				mat[j] ^= mat[i];
				inv[j] ^= inv[i];
			}
		}
	}

	/* Now the matrix is upper triangular.  Start at the top and multiply down */

	for (i = rows-1; i >= 0; i--) {
		for (j = 0; j < i; j++) {
			if (mat[j] & (one << i)) {
				inv[j] ^= inv[i];
			}
		}
	}

}

uint64_t galois_w64_bitmatrix_inverse (uint64_t x) {
	uint64_t mat[64];
	uint64_t inv[64];
	uint64_t ind, i, j, k, prod;
	uint64_t pp, one;

	one = 1;
	pp = 0;
	pp = 0x1b;

	for (i = 0; i < 64; i++) {
		mat[i] = x;

		if (x & (one << 63)) {
			x = x << 1;
			x = x ^ pp;
		} else {
			x = x << 1;
		}
	}

	galois_w64_invert_binary_matrix(mat, inv);

	return inv[0];
}

uint64_t galois_w64_shift_multiply(uint64_t x, uint64_t y) {
	uint64_t one, pp;
	uint64_t product[2];
	int i;

	product[0] = 0;
	product[1] = 0;
	one = 1;
	pp = 0x1b;
	for (i = 0; i < 64; i++) {
		if ((x & (one << i))) {
			product[0] ^= rshift(y, (64-i));
			product[1] ^= lshift(y, i);
			
			//product[0] ^= (y >> (64-i));
			//product[1] ^= (y << i);
		}
	}
	for (i = 127; i >= 64; i--) {
		//if ((product[0] & (one << (i-64)))) {
		if ((product[0] & lshift(one, i-64))) {

			product[0] ^= rshift(pp, (128-i));
			product[1] ^= lshift(pp, i-64);

			//product[0] ^= (pp >> (128-i));
			//product[1] ^= (pp << (i-64));
		}
	}
	return product[1];
}

int galois_w64_create_split8_tables () {
	uint64_t p1, p2, i, j, p1elt, p2elt, index, ishift, jshift, *table;
	uint64_t a, b;
	
	if (galois_w64_split8_tables[0] != NULL) return 0;

	for (i = 0; i < 15; i++) {
		galois_w64_split8_tables[i] = (uint64_t *) malloc(sizeof(uint64_t) * (1 << 8) * (1 << 8));
		if (galois_w64_split8_tables[i] == NULL) {
			for (i--; i >= 0; i--) free(galois_w64_split8_tables[i]);
			return -1;
		}
	}

	for (i = 0; i < 8; i+=7) {
		ishift = i * 8;
		for (j = (i == 0)? 0 : 1; j < 8; j++) {
			jshift = j * 8;
			table = galois_w64_split8_tables[i+j];
			index = 0;
			for (p1 = 0; p1 < (1 << 8); p1++) {
				p1elt = (p1 << ishift);
				for (p2 = 0; p2 < (1 << 8); p2++) {
					p2elt = (p2 << jshift);
					table[index] = galois_w64_bitmatrix_multiply(p1elt, p2elt);
					index++;
				}
			} 
		}
	}

	return 0;
}

uint64_t galois_w64_split8_table_multiply (uint64_t x, uint64_t y) {
	uint64_t a, b, accumulator, i8, j8;
	uint64_t i, j;
	uint64_t index_mod;

	accumulator = 0;
	index_mod = 255;

	i8 = 0;
	for (i = 0; i < 8; i++) {
		a = (((x >> i8) & index_mod) << 8);
		j8 = 0;
		for (j = 0; j < 8; j++) {
			b = ((y >> j8) & index_mod);
			accumulator ^= galois_w64_split8_tables[i+j][a|b];
			j8 += 8;
		}
		i8 += 8;
	}
	return accumulator;

}

uint64_t galois_w64_split8_table_divide (uint64_t x, uint64_t y) {
	uint64_t inverse;
	if (y == 0) return -1;
	if (x == 0) return 0;

	inverse = galois_w64_bitmatrix_inverse(y);
	return galois_w64_split8_table_multiply(x, inverse);
}

uint64_t galois_w64_composite_k2_multiply(uint64_t x, uint64_t y) {
	if (x == 0 || y == 0) return 0;
	int l;
	uint64_t s;
	uint64_t product, mask0;
	uint64_t p2, p3;
	uint64_t one;

	if (galois_w32_split8_tables[0] == NULL) {
		if (galois_w32_create_split8_tables() < 0) {
			fprintf(stderr, "galois_w64_composite_k2_multiply: Error creating tables for w=32\n");
			exit(1);
		}
	}

	l = 32;
	s = 2;
	one = 1;
	mask0 = (one << l) - 1;

	p2 = galois_w32_split8_table_multiply((x >> l), (y >> l));
	p3 = p2;
	p2 = galois_w32_split8_table_multiply(s, p2);
	
	product = galois_w32_split8_table_multiply((x >> l), (y & mask0));
	product ^= galois_w32_split8_table_multiply((x & mask0), (y >> l));
	product ^= p2;
	product <<= l;
	product ^= p3;
	product ^= galois_w32_split8_table_multiply((x & mask0), (y & mask0));
	
	return product;
	
}

uint64_t galois_w64_composite_k2_divide(uint64_t x, uint64_t y) {
	uint64_t inverse;
	if (y == 0) return -1;
	if (x == 0) return 0;

	inverse = galois_w64_composite_k2_inverse (y);
	return galois_w64_composite_k2_multiply (x, inverse);
	
}

uint64_t galois_w64_composite_k2_inverse(uint64_t x) {
	if (x == 0) return -1;
	int l;
	uint64_t s;
	uint64_t c;
	uint64_t mask0, inv, one;
	uint64_t a1i, a0i;
	int sum_log;

	if (galois_w32_split8_tables[0] == NULL) {
		if (galois_w32_create_split8_tables() < 0) {
			fprintf(stderr, "galois_w64_composite_k2_inverse: Error creating tables for w=32\n");
			exit(1);
		}
	}

	one = 1;
	l = 32;
	s = 2;
	mask0 = (one << l) - 1;

	if (((x >> l) & mask0) == 0) {
		inv = galois_bitmatrix_divide(1, (x & mask0), l);
		
	} else if ((x & mask0) == 0) {
		a1i = galois_bitmatrix_divide(1, ((x >> l) & mask0), l);
		inv = a1i;
		inv <<= l;
		inv ^= galois_w32_split8_table_multiply(a1i, s);
		
	} else {
		a1i = galois_bitmatrix_divide(1, ((x>>l) & mask0), l);
		a0i = galois_bitmatrix_divide(1, (x &mask0), l);
		
		c = galois_w32_split8_table_multiply(a0i, ((x >> l) & mask0));
		c^= galois_w32_split8_table_multiply((x & mask0), a1i);
		c^= (s & mask0);
		c = galois_bitmatrix_divide(1, (c & mask0), l);
		c = galois_w32_split8_table_multiply((c & mask0), ((x >> l) & mask0));
		c = galois_w32_split8_table_multiply(a0i, (c & mask0));
		
		inv = galois_w32_split8_table_multiply(a1i, (c & mask0));
		inv <<= l;
		inv^= galois_w32_split8_table_multiply(a0i, c^1);

	}
	return inv;
	
}

uint64_t galois_w64_recursive_composite_k2_multiply(uint64_t x, uint64_t y) {
	if (x == 0 || y == 0) return 0;
	uint64_t l;
	uint64_t s;
	uint64_t product, mask0;
	uint64_t p2, p3;
	uint64_t one;

	one = 1;
	l = 32;
	s = 2;

	mask0 = (one << l) - 1;

	p2 = galois_w32_composite_k2_multiply(x >> l, y >> l);
	p3 = p2;
	p2 = galois_w32_composite_k2_multiply(s, p2);

	product = galois_w32_composite_k2_multiply(x >> l, y & mask0);
	product ^= galois_w32_composite_k2_multiply(x & mask0, y >> l);
	product ^= p2;
	product <<= l;
	product ^= p3;
	product ^= galois_w32_composite_k2_multiply(x & mask0, y & mask0);

	return product;
	
}

uint64_t galois_w64_recursive_composite_k2_divide(uint64_t x, uint64_t y) {
	uint64_t inverse;
	if (y == 0) return -1;
	if (x == 0) return 0;

	inverse = galois_w64_recursive_composite_k2_inverse (y);
	return galois_w64_recursive_composite_k2_multiply (x, inverse);
	
}

uint64_t galois_w64_recursive_composite_k2_inverse(uint64_t x) {
	if (x == 0) return 0;
	uint64_t l;
	uint64_t s;
	uint64_t c;
	uint64_t mask0, inv;
	uint64_t a1i, a0i;
	uint64_t one;

	one = 1;
	l = 32;
	s = 2;
	mask0 = (one << l) - 1;

	//if (((x >> l) & mask0) == 0) {
	if ((rshift(x, l) & mask0) == 0) {
		inv = galois_w32_composite_k2_divide(1, x & mask0);
	} else if ((x & mask0) == 0) {
		//a1i = galois_w32_composite_k2_divide(1, (x >> l) & mask0);
		a1i = galois_w32_composite_k2_divide(1, rshift(x,l) & mask0);
		inv = a1i;
		inv <<= l;
		inv ^= galois_w32_composite_k2_multiply(a1i, s);
	} else {
		//a1i = galois_w32_composite_k2_divide(1, (x >> l) & mask0);
		a1i = galois_w32_composite_k2_divide(1, rshift(x,l) & mask0);
		a0i = galois_w32_composite_k2_divide(1, x & mask0);

		//c = galois_w32_composite_k2_multiply((x >> l) & mask0, a0i);
		c = galois_w32_composite_k2_multiply(rshift(x,l) & mask0, a0i);
		c ^= galois_w32_composite_k2_multiply(x & mask0, a1i);
		c ^= s & mask0;
		c = galois_w32_composite_k2_divide(1, c & mask0);
		//c = galois_w32_composite_k2_multiply(c & mask0, (x >> l) & mask0);
		c = galois_w32_composite_k2_multiply(c & mask0, rshift(x,l) & mask0);
		c = galois_w32_composite_k2_multiply(c & mask0, a0i);

		inv = galois_w32_composite_k2_multiply(a1i, c & mask0);
		inv <<= l;
		inv ^= galois_w32_composite_k2_multiply(c ^ 1, a0i);
	}
	return inv;

}

/* w == 128 */

uint64_t *galois_w128_bitmatrix_multiply (uint64_t *x, uint64_t *y) {
	int i, k;
	uint64_t *prod;
	uint64_t scratch[128][2];
	uint64_t *s, one, pp, j, wshift;
	uint64_t savey[2];

	prod = (uint64_t *) malloc (sizeof(uint64_t) * 2);
	one = 1;
	pp = 0x87;
	savey[0] = y[0]; savey[1] = y[1];
	wshift = lshift(one, 63);

	for (i = 0; i < 128; i++) {
		scratch[i][0] = savey[0];
		scratch[i][1] = savey[1];

		if (savey[0] & wshift) {
			savey[0] = savey[0] << 1;
			if ((savey[1] & wshift) != 0) savey[0] = savey[0] ^ 1;
			savey[1] = savey[1] << 1;
			savey[1] = savey[1] ^ pp;
		} else {
			savey[0] = savey[0] << 1;
			if ((savey[1] & wshift) != 0) savey[0] = savey[0] ^ 1;
			savey[1] = savey[1] << 1;
		}
	}
	
	prod[0] = 0; prod[1] = 0;
	/* Iterate over the first 64 bits of x */
	for (i = 0; i < 64; i++) {
		if (x[1] & (one << i)) {
			prod[1] = prod[1] ^ scratch[i][1];
			prod[0] = prod[0] ^ scratch[i][0];
		}
		if (x[0] & (one << i)) {
			prod[1] = prod[1] ^ scratch[i+64][1];
			prod[0] = prod[0] ^ scratch[i+64][0];
		}
	}

	return prod;
}

uint64_t *galois_w128_bitmatrix_divide (uint64_t *x, uint64_t *y) {
	uint64_t *inverse;
	uint64_t *retval;
	retval = (uint64_t *) malloc(sizeof(uint64_t) * 2);
	retval[0] = 0;
	if (y[0] == 0 && y[1] == 0) { retval[1] = -1; return retval; }
	if (x[0] == 0 && x[1] == 0) { retval[1] = 0;  return retval; }
	free(retval);

	inverse = galois_w128_bitmatrix_inverse(y);
	retval = galois_w128_bitmatrix_multiply(x, inverse);
	free(inverse);
	return retval;
}

void galois_w128_invert_binary_matrix (uint64_t **mat, uint64_t **inv) {
	uint64_t cols, rows;
	int i, j, k;
	uint64_t tmp[2];
	uint64_t one;

	one = 1;
	cols = 128;
	rows = 128;

	for (i = 0; i < rows; i++) {
		inv[i][0] = 0;
		inv[i][1] = 0;
		if (i < 64) {
			inv[i][1] = (one << i);
		} else {
			inv[i][0] = (one << (i - 64));
		}
	}

	/* First -- convert into upper triangular */

	for (i = 0; i < cols; i++) {

		/* Swap rows if we ave a zero i,i element.  If we can't swap, then the
			 matrix was not invertible */

		if (i < 64) {
			if ((mat[i][1] & (one << i)) == 0) {
				for (j = i+1; j < rows && (mat[j][1] & (one << i)) == 0; j++);
				if (j == rows) {
					fprintf(stderr, "galois_invert_matrix: Matrix not invertible!! %d\n", i);
					exit(1);
				}
				tmp[0] = mat[i][0]; mat[i][0] = mat[j][0]; mat[j][0] = tmp[0];
				tmp[1] = mat[i][1]; mat[i][1] = mat[j][1]; mat[j][1] = tmp[1];
				tmp[0] = inv[i][0]; inv[i][0] = inv[j][0]; inv[j][0] = tmp[0];
				tmp[1] = inv[i][1]; inv[i][1] = inv[j][1]; inv[j][1] = tmp[1];
			}
		} else {
			if ((mat[i][0] & (one << (i-64))) == 0) {
				for (j = i+1; j < rows && (mat[j][0] & (one << (i-64))) == 0; j++);
				if (j == rows) {
					fprintf(stderr, "galois_invert_matrix: Matrix not invertible!!\n");
					exit(1);
				}
				tmp[0] = mat[i][0]; mat[i][0] = mat[j][0]; mat[j][0] = tmp[0];
				tmp[1] = mat[i][1]; mat[i][1] = mat[j][1]; mat[j][1] = tmp[1];
				tmp[0] = inv[i][0]; inv[i][0] = inv[j][0]; inv[j][0] = tmp[0];
				tmp[1] = inv[i][1]; inv[i][1] = inv[j][1]; inv[j][1] = tmp[1];
			}
		}

		/* Now for each j>i, add A_ji*Ai to Aj */
		for (j = i+1; j != rows; j++) {
			if (i < 64 && (mat[j][1] & (one << i)) != 0) {
				mat[j][1] ^= mat[i][1];
				mat[j][0] ^= mat[i][0];
				inv[j][1] ^= inv[i][1];
				inv[j][0] ^= inv[i][0];
			} else if (i >= 64 && (mat[j][0] & (one << (i-64))) != 0) {
				mat[j][1] ^= mat[i][1];
				mat[j][0] ^= mat[i][0];
				inv[j][1] ^= inv[i][1];
				inv[j][0] ^= inv[i][0];
			}
		}
	}

	/* Now the matrix is upper triangular.  Start at the top and multiply down */

	for (i = rows-1; i >= 0; i--) {
		for (j = 0; j < i; j++) {
			if (i < 64 && (mat[j][1] & (one << i)) != 0) {
				inv[j][1] ^= inv[i][1];
				inv[j][0] ^= inv[i][0];
			} else if (i >= 64 && (mat[j][0] & (one << (i-64))) != 0) {
				inv[j][1] ^= inv[i][1];
				inv[j][0] ^= inv[i][0];
			}
		}
	}

}

uint64_t *galois_w128_bitmatrix_inverse (uint64_t *x) {
	uint64_t *mat[128];
	uint64_t *inv[128];
	uint64_t ind, i, j, k, prod;
	uint64_t pp, one;
	uint64_t savex[2];
	uint64_t *retval;

	savex[0] = x[0];
	savex[1] = x[1];

	one = 1;
	pp = 0;
	pp = 0x87;

	for (i = 0; i < 128; i++) {
		mat[i] = (uint64_t *) malloc(sizeof(uint64_t) * 2);
		inv[i] = (uint64_t *) malloc(sizeof(uint64_t) * 2);
	}

	for (i = 0; i < 128; i++) {
		mat[i][0] = savex[0];
		mat[i][1] = savex[1];

		if (savex[0] & (one << 63)) {
			savex[0] = savex[0] << 1;
			if (savex[1] & (one << 63)) savex[0] ^= 1;
			savex[1] = savex[1] << 1;
			savex[1] = savex[1] ^ pp;
		} else {
			savex[0] = savex[0] << 1;
			if (savex[1] & (one << 63)) savex[0] ^= 1;
			savex[1] = savex[1] << 1;
		}
	}

	galois_w128_invert_binary_matrix (mat, inv);
	retval = (uint64_t *) malloc (sizeof(uint64_t) * 2);
	retval[0] = inv[0][0];
	retval[1] = inv[0][1];
	for (i = 0; i < 128; i++) {
		free(mat[i]);
		free(inv[i]);
	}

	return retval;
}

int galois_w128_create_split8_tables() {
	uint64_t p1, p2, i, j, p1elt[2], p2elt[2], index, ishift, jshift, *table, *tmp;
	uint64_t a, b;

	if (galois_w128_split8_tables[0] != NULL) return 0;
	for (i = 0; i < 31; i++) {
		galois_w128_split8_tables[i] = (uint64_t *) malloc(sizeof(uint64_t) * (1 << 8) * (1 << 8) * 2);
		if (galois_w128_split8_tables[i] == NULL) {
			for (i--; i >= 0; i--) free(galois_w128_split8_tables[i]);
			return -1;
		}
	}

	for (i = 0; i < 16; i+=15) {
		ishift = i * 8;
		for (j = (i == 0)? 0 : 1; j < 16; j++) {
			jshift = j * 8;
			table = galois_w128_split8_tables[i+j];
			index = 0;
			for (p1 = 0; p1 < (1 << 8); p1++) {
				if (i != 0) {
					p1elt[0] = lshift(p1, ishift-64);
					p1elt[1] = 0;
				} else {
					p1elt[0] = 0;
					p1elt[1] = p1;
				}

				for (p2 = 0; p2 < (1 << 8); p2++) {
					if (jshift < 64) {
						p2elt[0] = 0;
						p2elt[1] = lshift(p2, jshift);
					} else {
						p2elt[0] = lshift(p2, jshift-64);
						p2elt[1] = 0;
					}

					tmp = galois_w128_bitmatrix_multiply(p1elt, p2elt);
					table[index] = tmp[0];
					table[index+1] = tmp[1];
					index+=2;
					free(tmp);
				}
			} 
		}
	}

	return 0;
}

uint64_t *galois_w128_split8_table_multiply(uint64_t *x, uint64_t *y) {
	int i, j, a, b, i8, j8;
	uint64_t *accumulator;

	accumulator = (uint64_t *) malloc(sizeof(uint64_t) * 2);
	accumulator[0] = 0;
	accumulator[1] = 0;

	i8 = 0;
	for (i = 0; i < 16; i++) {
		if (i8 >= 64) a = (((rshift(x[0], i8-64)) & 255) << 9);
		else a = (((rshift(x[1], i8)) & 255) << 9);
		j8 = 0;
		for (j = 0; j < 16; j++) {
			if (j8 >= 64) b = (((rshift(y[0], j8-64)) & 255) << 1);
			else b = (((rshift(y[1], j8)) & 255) << 1);
			accumulator[0] ^= galois_w128_split8_tables[i+j][a|b];
			accumulator[1] ^= galois_w128_split8_tables[i+j][(a|b) + 1];
			j8 += 8;
		}
		i8 += 8;
	}
	return accumulator;

}

uint64_t *galois_w128_shift_multiply(uint64_t *x, uint64_t *y) {
	uint64_t one, pp;
	uint64_t product[4];
	uint64_t *retval;
	int i;

	product[0] = 0;
	product[1] = 0;
	product[2] = 0;
	product[3] = 0;
	
	one = 1;
	pp = 0x87;
	for (i = 0; i < 128; i++) {
		if ((i<64) && (x[1] & lshift(one, i))) {
			product[1] ^= rshift(y[0], 64-i);
			product[2] ^= lshift(y[0], i);
			product[2] ^= rshift(y[1], 64-i);
			product[3] ^= lshift(y[1], i);
			/*
			product[1] ^= (y[0] >> (64-i));
			product[2] ^= (y[0] << i);
			product[2] ^= (y[1] >> (64-i));
			product[3] ^= (y[1] << i);
			*/
		} else if (i >= 64 && (x[0] & lshift(one, i-64))) {
			product[0] ^= rshift(y[0], 128-i);
			product[1] ^= lshift(y[0], i-64);
			product[1] ^= rshift(y[1], 128-i);
			product[2] ^= lshift(y[1], i-64);
			/*
			product[0] ^= (y[0] >> (128 -i));
			product[1] ^= (y[0] << (i - 64));
			product[1] ^= (y[1] >> (128 - i));
			product[2] ^= (y[1] << (i-64));
			*/
		}
	}
	for (i = 255; i >= 128; i--) {
		if ((i >= 192) && (product[0] & (lshift(one, i-192)))) {
			product[1] ^= rshift(pp, 64-(i-192));
			product[2] ^= lshift(pp, i-192);
			/*
			product[1] ^= (pp >> (64-(i-192)));
			product[2] ^= (pp << (i-192));
			*/
		} else if ((i >= 128) && (product[1] & (lshift(one, i-128)))) {
			product[2] ^= rshift(pp, 64-(i-128));
			product[3] ^= lshift(pp, i-128);
			/*
			product[2] ^= (pp >> (64-(i-128)));
			product[3] ^= (pp << (i-128));
			*/
		}
	}
	retval = (uint64_t *) malloc (sizeof(uint64_t) * 2);
	retval[0] = product[2];
	retval[1] = product[3];
	return retval;
}

uint64_t *galois_w128_composite_k2_multiply(uint64_t *x, uint64_t *y) {
	int l;
	uint64_t s;
	uint64_t *product;
	uint64_t p2, p3;
	uint64_t one;
	
	if (x[0] == 0 && x[1] == 0) {
		product = talloc(uint64_t, 2);
		product[1] = 0;
		product[0] = 0;
		return product;
	}
	
	if (y[1] == 0 && y[1] == 0) {
		product = talloc(uint64_t, 2);
		product[1] = 0;
		product[0] = 0;
		return product;
	}
	

	if (galois_w64_split8_tables[0] == NULL) {
		if (galois_w64_create_split8_tables() < 0) {
			fprintf(stderr, "galois_w128_composite_k2_multiply: Error creating tables for w=64\n");
			exit(1);
		}
	}

	product = talloc(uint64_t, 2);

	l = 64;
	s = 2;
	one = 1;

	p2 = galois_w64_split8_table_multiply(x[0], y[0]);
	p3 = p2;
	p2 = galois_w64_split8_table_multiply(s, p2);
	
	product[1]  = galois_w64_split8_table_multiply(x[0], y[1]);
	product[1] ^= galois_w64_split8_table_multiply(x[1], y[0]);
	product[1] ^= p2;
	product[0]  = product[1];
	product[1]  = p3;
	product[1] ^= galois_w64_split8_table_multiply(x[1], y[1]);
	
	return product;

}

uint64_t *galois_w128_composite_k2_divide(uint64_t *x, uint64_t *y) {
	uint64_t *inverse, *retval;
	if (y == 0) return NULL;
	if (x == 0) {
		retval  = talloc (uint64_t, 2);
		retval[0] = 0;
		retval[1] = 0;
		return retval;
	}

	inverse = galois_w128_composite_k2_inverse (y);
	retval = galois_w128_composite_k2_multiply (x, inverse);
	free (inverse);
	return retval;

}

uint64_t *galois_w128_composite_k2_inverse(uint64_t *x) {
	if (x[0] == 0 && x[1] == 0) {
		return NULL;
	}
	int l;
	uint64_t s;
	uint64_t c;
	uint64_t *inv, one;
	uint64_t a1i, a0i;

	if (galois_w64_split8_tables[0] == NULL) {
		if (galois_w64_create_split8_tables() < 0) {
			fprintf(stderr, "galois_w128_composite_k2_inverse: Error creating tables for w=64\n");
			exit(1);
		}
	}

	inv = talloc (uint64_t, 2);
	inv[0] = 0;
	inv[1] = 0;
	one = 1;
	l = 64;
	s = 2;

	if (x[0] == 0) {
		inv[1] = galois_w64_bitmatrix_divide(1, x[1]);
		
	} else if (x[1] == 0) {
		a1i    = galois_w64_bitmatrix_divide(1, x[0]);
		inv[0] = a1i;
		inv[1] = galois_w64_split8_table_multiply(a1i, s);
		
	} else {
		a1i = galois_w64_bitmatrix_divide(1, x[0]);
		a0i = galois_w64_bitmatrix_divide(1, x[1]);
		
		c = galois_w64_split8_table_multiply(a0i, x[0]);
		c^= galois_w64_split8_table_multiply(x[1], a1i);
		c^= s;
		c = galois_w64_bitmatrix_divide(1, c);
		c = galois_w64_split8_table_multiply(c, x[0]);
		c = galois_w64_split8_table_multiply(a0i, c);
		
		inv[0] = galois_w64_split8_table_multiply(a1i, c);
		inv[1] = galois_w64_split8_table_multiply(a0i, c^1);

	}
	return inv;

}

uint64_t *galois_w128_recursive_composite_k2_multiply(uint64_t *x, uint64_t *y) {
	uint64_t s;
	uint64_t *product;
	uint64_t p2, p3;
	
	if (x[0] == 0 && x[1] == 0) {
		product = talloc(uint64_t, 2);
		product[1] = 0;
		product[0] = 0;
		return product;
	}
	
	if (y[1] == 0 && y[1] == 0) {
		product = talloc(uint64_t, 2);
		product[1] = 0;
		product[0] = 0;
		return product;
	}

	s = 2;
	product = talloc(uint64_t, 2);

	p2 = galois_w64_recursive_composite_k2_multiply(x[0], y[0]);
	p3 = p2;
	p2 = galois_w64_recursive_composite_k2_multiply(s, p2);

	product[0] = galois_w64_recursive_composite_k2_multiply(x[0], y[1]);
	product[0] ^= galois_w64_recursive_composite_k2_multiply(x[1], y[0]);
	product[0] ^= p2;

	product[1]  = p3;
	product[1] ^= galois_w64_recursive_composite_k2_multiply(x[1], y[1]);

	return product;

}

uint64_t *galois_w128_recursive_composite_k2_divide(uint64_t *x, uint64_t *y) {
	uint64_t *inverse, *retval;
	if (y[1] == 0 && y[0] == 0) return NULL;
	if (x[0] == 0 && x[1] == 0) {
		retval = talloc(uint64_t, 2);
		retval[0] = 0;
		retval[1] = 0;
		return retval;
	}

	inverse = galois_w128_recursive_composite_k2_inverse (y);
	retval = galois_w128_recursive_composite_k2_multiply (x, inverse);
	free(inverse);
	return retval;

}

uint64_t *galois_w128_recursive_composite_k2_inverse(uint64_t *x) {
	uint64_t s;
	uint64_t c;
	uint64_t *inv;
	uint64_t a1i, a0i;

	if (x[0] == 0 && x[1] == 0) return NULL;

	s = 2;

	inv = talloc(uint64_t, 2);
	inv[1] = 0;
	inv[0] = 0;

	if (x[0] == 0) {
		inv[1] = galois_w64_recursive_composite_k2_divide(1, x[1]);
	} else if (x[1] == 0) {
		a1i = galois_w64_recursive_composite_k2_divide(1, x[0]);
		inv[0] = a1i;
		inv[1] = galois_w64_recursive_composite_k2_multiply(a1i, s);
	} else {
		a1i = galois_w64_recursive_composite_k2_divide(1, x[0]);
		a0i = galois_w64_recursive_composite_k2_divide(1, x[1]);

		c = galois_w64_recursive_composite_k2_multiply(x[0], a0i);
		c ^= galois_w64_recursive_composite_k2_multiply(x[1], a1i);
		c ^= s;
		c = galois_w64_recursive_composite_k2_divide(1, c);
		c = galois_w64_recursive_composite_k2_multiply(c, x[0]);
		c = galois_w64_recursive_composite_k2_multiply(c, a0i);

		inv[0] = galois_w64_recursive_composite_k2_multiply(a1i, c);
		inv[1] = galois_w64_recursive_composite_k2_multiply(c ^ 1, a0i);
	}
	return inv;

}

/* access functions */

void *galois_get_mult_table(int w) {
	if (w <= 8) {
		if (galois_w08_mult_tables[w] == NULL) {
			if (galois_w08_create_mult_tables(w) < 0) {
				return NULL;
			}
		}
		return (void *) galois_w08_mult_tables[w];

	} else if (w <= 16) {
		if (galois_w16_mult_tables[w] == NULL) {
			if (galois_w16_create_mult_tables(w) < 0) {
				return NULL;
			}
		}
		return (void *) galois_w16_mult_tables[w];

	} else {
		fprintf(stderr, "Galois_get_mult_table - no implementation for w=%d\n", w);
		exit(1);
	}
}

void *galois_get_div_table(int w)  {
	if (w <= 8) {
		if (galois_w08_div_tables[w] == NULL) {
			if (galois_w08_create_mult_tables(w) < 0) {
				return NULL;
			}
		}
		return (void *) galois_w08_div_tables[w];

	} else if (w <= 16) {
		if (galois_w16_div_tables[w] == NULL) {
			if (galois_w16_create_mult_tables(w) < 0) {
				return NULL;
			}
		}
		return galois_w16_div_tables[w];

	} else {
		fprintf(stderr, "Galois_get_div_table - no implementation for w=%d\n", w);
		exit(1);
	}
}

void *galois_get_log_table(int w) {
	if (w <= 8) {
		if (galois_w08_log_tables[w] == NULL) {
			if (galois_w08_create_log_tables(w) < 0) {
				return NULL;
			}
		}
		return (void *) galois_w08_log_tables[w];

	} else if (w <= 16) {
		if (galois_w16_log_tables[w] == NULL) {
			if (galois_w16_create_log_tables(w) < 0) {
				return NULL;
			}
		}
		return (void *) galois_w16_log_tables[w];

	} else if (w <= 32) {
		if (galois_w32_log_tables[w] == NULL) {
			if (galois_w32_create_log_tables(w) < 0) {
				return NULL;
			}
		}
		return (void *) galois_w32_log_tables[w];

	} else {
		fprintf(stderr, "Galois_get_log_table - no implementation for w=%d\n", w);
		exit(1);
	}
}

void *galois_get_ilog_table(int w) {
	if (w <= 8) {
		if (galois_w08_ilog_tables[w] == NULL) {
			if (galois_w08_create_log_tables(w) < 0) {
				return NULL;
			}
		}
		return (void *) galois_w08_ilog_tables[w];

	} else if (w <= 16) {
		if (galois_w16_ilog_tables[w] == NULL) {
			if (galois_w16_create_log_tables(w) < 0) {
				return NULL;
			}
		}
		return (void *) galois_w16_ilog_tables[w];

	} else if (w <= 32) {
		if (galois_w32_ilog_tables[w] == NULL) {
			if (galois_w32_create_log_tables(w) < 0) {
				return NULL;
			}
		}
		return (void *) galois_w32_ilog_tables[w];

	} else {
		fprintf(stderr, "Galois_get_ilog_table - no implementation for w=%d\n", w);
		exit(1);
	}
}

uint32_t **galois_get_split8_tables() { return galois_w32_split8_tables; }

/* Region multiply functions, organized by w */
/*  w=4 implementations:
	full mult tables
	lazy multtable
	shift
	bitmatrix
	log/antilog
*/

int galois_w04_create_shifted_mult_table () {
	int j, x, y, logx;

	if (galois_w08_mult_tables[4] == NULL) {
		if (galois_w08_create_mult_tables(4) < 0) {
			return -1;
		}
	}

	if (galois_w04_shifted_mult_table != NULL) return 0;
	galois_w04_shifted_mult_table = talloc (unsigned char, nw[4] * nw[4]);
	if (galois_w04_shifted_mult_table == NULL) return -1;

	/* Set mult/div tables for x = 0 */
	for (j = 0; j < (nw[4] * nw[4]); j++) {
		galois_w04_shifted_mult_table[j] = (galois_w08_mult_tables[4][j] << 4);
	}

	return 0;
}

void galois_w04_multtable_region_multiply(unsigned char *region,
		int multby,
		int nbytes,
		unsigned char *r2,  
		int add)
{
	unsigned char *ur1, *ur2, *cp;
	unsigned char prod;
	int i, srow, j;
	uint64_t l, *lp2;
	unsigned char *lp;
	int sol;

	ur1 = (unsigned char *) region;
	ur2 = (r2 == NULL) ? ur1 : (unsigned char *) r2;

	if (galois_w08_mult_tables[4] == NULL) {
		if (galois_w08_create_mult_tables(4) < 0) {
			fprintf(stderr, "galois_w04_multtable_region_multiply -- couldn't make multiplication tables\n");
			exit(1);
		}
	}
	if (galois_w04_shifted_mult_table == NULL) {
		if (galois_w04_create_shifted_mult_table() < 0) {
			fprintf(stderr, "galois_w04_multtable_region_multiply -- couldn't make shifted multiplication tables\n");
			exit(1);
		}
	}
	srow = multby * nw[4];
	if (r2 == NULL || !add) {
		for (i = 0; i < nbytes; i++) {
			prod = galois_w08_mult_tables[4][srow+ur1[i]];
			prod ^= galois_w04_shifted_mult_table[srow+ur1[i]];
			ur2[i] = prod;
		}
	} else {
		sol = sizeof(uint64_t);
		lp2 = &l;
		lp = (unsigned char *) lp2;
		for (i = 0; i < nbytes; i += sol) {
			cp = ur2+i;
			lp2 = (uint64_t *) cp;
			for (j = 0; j < sol; j++) {
				prod = galois_w08_mult_tables[4][srow+ur1[i+j]];
				prod ^= galois_w04_shifted_mult_table[srow+ur1[i+j]];
				lp[j] = prod;
			}
			*lp2 = (*lp2) ^ l;
		}
	}
	return;
	
}

unsigned char *galois_w04_create_lazytable_row (int row_index) {
	unsigned char *row_elts, tmp;
	int i;
	int logx;

	if (galois_w08_log_tables[4] == NULL) {
		if (galois_w08_create_log_tables(4) < 0) { return NULL; }
	}

	row_elts = (unsigned char *) malloc(sizeof(unsigned char) * nw[4] * nw[4]);

	logx = galois_w08_log_tables[4][row_index];
	row_elts[0] = 0;
	for (i = 1; i < nw[4]; i++) {
		row_elts[i] = (unsigned char) galois_w08_ilog_tables[4][logx + galois_w08_log_tables[4][i]];
		tmp = (unsigned char) galois_w08_ilog_tables[4][logx + galois_w08_log_tables[4][i]];
		row_elts[i+nw[4]] = (tmp << 4);
	}
	return row_elts;
}

void galois_w04_lazytable_region_multiply(unsigned char *region,
		int multby,
		int nbytes,
		unsigned char *r2,
		int add)
{
	unsigned char *table_row, *shifted_row;
	unsigned char *ur1, *ur2, *cp;
	unsigned char prod;
	int i, srow, j;
	uint64_t l, *lp2, *lptop;
	unsigned char *lp;
	int sol;

	ur1 = (unsigned char *) region;
	ur2 = (r2 == NULL) ? ur1 : (unsigned char *) r2;

	if (multby == 0) {
		bzero(ur2, nbytes);
		return;
	}

	table_row = galois_w04_create_lazytable_row(multby);
	if (table_row == NULL) {
		fprintf(stderr, "galois_08_lazytable_region_multiply -- couldn't make multiplication tables\n");
		exit(1);
	}
	shifted_row = table_row + nw[4];

	if (r2 == NULL || !add) {
		for (i = 0; i < nbytes; i++) {
			prod = table_row[ur1[i]];
			prod ^= shifted_row[ur1[i]];
			ur2[i] = prod;
		}
	} else {
		sol = sizeof(uint64_t);
		lp2 = &l;
		lp = (unsigned char *) lp2;
		for (i = 0; i < nbytes; i += sol) {
			cp = ur2+i;
			lp2 = (uint64_t *) cp;
			for (j = 0; j < sol; j++) {
				prod = table_row[ur1[i+j]];
				prod^= shifted_row[ur1[i+j]];
				lp[j] = prod;
			}
			*lp2 = (*lp2) ^ l;
		}
	}
	free (table_row);
	return;

}

void galois_w04_shift_region_multiply(unsigned char *region,
		int multby,
		int nbytes,
		unsigned char *r2, 
		int add)
{
	unsigned char *ur1, *ur2, *ur2top;
	uint16_t product, pp, pps, n1;
	int i, j;
	int num_ones, scratch[4], index, top;

	ur1 = (unsigned char *) region;
	ur2 = (r2 == NULL) ? ur1 : (unsigned char *) r2;
	ur2top = ur2 + nbytes;

	if (multby == 0) {
		if (!add) {
			while (ur2 < ur2top) {
				*ur2 = 0;
				ur2++;
			}
		}
		return;
	}

	pp = prim_poly[4];
	pps = (prim_poly[4] << 8);
	top = 0;
	index = 0;
	num_ones = 0;
	for (i = 0; i < 4; i++) {
		if (multby & (1 << i)) {
			scratch[index] = i;
			index++;
			num_ones++;
			top = i;
		}
	}
	if (!add) {
		for (j = 0; j < nbytes; j++) {
			product = 0;
			n1  = (ur1[j] & 0xf0) << 4;
			n1 += (ur1[j] & 0x0f);
			for (i = 0; i < num_ones; i++) product ^= (n1 << scratch[i]);
			for (i = top+12; i >= 8; i--) if(product & (1 << i)) product ^= (pps << (i - 4));
			for (i = top+4; i>= 4; i--) if (product & (1 << i)) product ^= (pp << (i-4));
			i = (product & 0xf00) >> 4;
			i+= (product & 0xf);
			ur2[j] = i;
		}
	} else {
		for (j = 0; j < nbytes; j++) {
			product = 0;
			n1  = (ur1[j] & 0xf0) << 4;
			n1 += (ur1[j] & 0x0f);
			for (i = 0; i < num_ones; i++) product ^= (n1 << scratch[i]);
			for (i = top+12; i >= 8; i--) if(product & (1 << i)) product ^= (pps << (i - 4));
			for (i = top+4; i>= 4; i--) if (product & (1 << i)) product ^= (pp << (i-4));
			i = (product & 0xf00) >> 4;
			i+= (product & 0xf);
			ur2[j] ^= i;
		}
	}
	return;

}

void galois_w04_bitmatrix_region_multiply(unsigned char *region,
		int multby,
		int nbytes,
		unsigned char *r2, 
		int add)
{
	unsigned char *ur1, *ur2, *ur2top;
	unsigned char prod;
	uint16_t pp;
	int i, j;
	int scratch[4], shifted_scratch[4];

	ur1 = (unsigned char *) region;
	ur2 = (r2 == NULL) ? ur1 : (unsigned char *) r2;
	ur2top = ur2 + nbytes;

	if (multby == 0) {
		if (!add) {
			while (ur2 < ur2top) {
				*ur2 = 0;
				ur2++;
			}
		}
		return;
	}

	pp = prim_poly[4];

	for (i = 0; i < 4; i++) {
		scratch[i] = multby;
		shifted_scratch[i] = (multby << 4);
		multby <<= 1;
		if ((multby & nw[4])) {
			multby ^= pp;
		}
	}

	if (!add) {
		for (j = 0; j < nbytes; j++) {
			prod = 0;
			for (i = 0; i < 4; i++) {
				if ((ur1[j] & (1 << i))) prod ^= scratch[i];
			}
			for (; i < 8; i++) {
				if ((ur1[j] & (1 << i))) prod ^= shifted_scratch[i-4];
			}
			ur2[j] = prod;
		}
	} else {
		for (j = 0; j < nbytes; j++) {
			prod = 0;
			for (i = 0; i < 4; i++) {
				if ((ur1[j] & (1 << i))) prod ^= scratch[i];
			}
			for (; i < 8; i++) {
				if ((ur1[j] & (1 << i))) prod ^= shifted_scratch[i-4];
			}
			ur2[j] ^= prod;
		}
	}
	return;

}

void galois_w04_logtable_region_multiply(unsigned char *region,
		int multby,
		int nbytes,
		unsigned char *r2,
		int add)
{
	unsigned char *ur1, *ur2, *cp;
	int prod;
	int i, log1, j, log2;
	uint64_t l, *lp2, *lptop;
	unsigned char *lp;
	int sol;

	ur1 = (unsigned char *) region;
	ur2 = (r2 == NULL) ? ur1 : (unsigned char *) r2;

	if (multby == 0) {
		if (!add) {
			lp2 = (uint64_t *) ur2;
			ur2 += nbytes;
			lptop = (uint64_t *) ur2;
			while (lp2 < lptop) { *lp2 = 0; lp2++; }
		}
		return;
	}

	if (galois_w08_log_tables[4] == NULL) {
		if (galois_create_log_tables(4) < 0) {
			fprintf(stderr, "galois_w08_logtable_region_multiply -- couldn't make log tables\n");
			exit(1);
		}
	}
	log1 = galois_w08_log_tables[4][multby];

	if (r2 == NULL || !add) {
		for (i = 0; i < nbytes; i++) {
			if (ur1[i] == 0) {
				ur2[i] = 0;
			} else {
				log2 = galois_w08_log_tables[4][(ur1[i] & 0xf)] + log1;
				prod = (galois_w08_ilog_tables[4][log2] << 4);
				log2 = galois_w08_log_tables[4][(ur1[i] & 0xf0)] + log1;
				prod += (galois_w08_ilog_tables[4][log2]);
				ur2[i] = prod;
			}
		}
	} else {
		sol = sizeof(uint64_t);
		lp2 = &l;
		lp = (unsigned char *) lp2;
		for (i = 0; i < nbytes; i += sol) {
			cp = ur2+i;
			lp2 = (uint64_t *) cp;
			for (j = 0; j < sol; j++) {
				if (ur1[i+j] == 0) {
					lp[j] = 0;
				} else {
					log2 = galois_w08_log_tables[4][(ur1[i+j] & 0xf)] + log1;
					prod = (galois_w08_ilog_tables[4][log2] << 4);
					log2 = galois_w08_log_tables[4][(ur1[i+j] & 0xf0)] + log1;
					prod += (galois_w08_ilog_tables[4][log2]);
					lp[j] = prod;
				}
			}
			*lp2 = (*lp2) ^ l;
		}
	}
	return;
}

/*  w=8 implementations:
	full mult tables
	lazy multtable (only one row in memory at a time)
	shift
	bitmatrix
	log/antilog
	composite k=2 */

void galois_w08_multtable_region_multiply(unsigned char *region,      	/* Region to multiply */
		int multby,       											  	/* Number to multiply by */
		int nbytes,        											  	/* Number of bytes in region */
		unsigned char *r2,          								  	/* If r2 != NULL, products go here */
		int add) 
{
	unsigned char *ur1, *ur2, *cp;
	unsigned char prod;
	int i, srow, j;
	uint64_t l, *lp2;
	unsigned char *lp;
	int sol;

	ur1 = (unsigned char *) region;
	ur2 = (r2 == NULL) ? ur1 : (unsigned char *) r2;

	if (galois_w08_mult_tables[8] == NULL) {
		if (galois_w08_create_mult_tables(8) < 0) {
			fprintf(stderr, "galois_w08_multtable_region_multiply -- couldn't make multiplication tables\n");
			exit(1);
		}
	}
	srow = multby * nw[8];
	if (r2 == NULL || !add) {
		for (i = 0; i < nbytes; i++) {
			prod = galois_w08_mult_tables[8][srow+ur1[i]];
			ur2[i] = prod;
		}
	} else {
		sol = sizeof(uint64_t);
		lp2 = &l;
		lp = (unsigned char *) lp2;
		for (i = 0; i < nbytes; i += sol) {
			cp = ur2+i;
			lp2 = (uint64_t *) cp;
			for (j = 0; j < sol; j++) {
				prod = galois_w08_mult_tables[8][srow+ur1[i+j]];
				lp[j] = prod;
			}
			*lp2 = (*lp2) ^ l;
		}
	}
	return;
}

unsigned char *galois_w08_create_lazytable_row (int row_index) {
	unsigned char *row_elts;
	int i;
	int logx;

	if (galois_w08_log_tables[8] == NULL) {
		if (galois_w08_create_log_tables(8) < 0) { return NULL; }
	}

	row_elts = (unsigned char *) malloc(sizeof(unsigned char) * (1 << 8));

	logx = galois_w08_log_tables[8][row_index];
	row_elts[0] = 0;
	for (i = 1; i < (1 << 8); i++) {
		row_elts[i] = (unsigned char) galois_w08_ilog_tables[8][logx + galois_w08_log_tables[8][i]];
	}
	return row_elts;
}

void galois_w08_lazytable_region_multiply (unsigned char *region,
		int multby,
		int nbytes,
		unsigned char *r2,
		int add)
{
	unsigned char *table_row;
	unsigned char *ur1, *ur2, *cp;
	unsigned char prod;
	int i, srow, j;
	uint64_t l, *lp2, *lptop;
	unsigned char *lp;
	int sol;

	ur1 = (unsigned char *) region;
	ur2 = (r2 == NULL) ? ur1 : (unsigned char *) r2;

	if (multby == 0) {
		bzero(ur2, nbytes);
		return;
	}

	table_row = galois_w08_create_lazytable_row(multby);
	if (table_row == NULL) {
		fprintf(stderr, "galois_08_lazytable_region_multiply -- couldn't make multiplication tables\n");
		exit(1);
	}

	if (r2 == NULL || !add) {
		for (i = 0; i < nbytes; i++) {
			prod = table_row[ur1[i]];
			ur2[i] = prod;
		}
	} else {
		sol = sizeof(uint64_t);
		lp2 = &l;
		lp = (unsigned char *) lp2;
		for (i = 0; i < nbytes; i += sol) {
			cp = ur2+i;
			lp2 = (uint64_t *) cp;
			for (j = 0; j < sol; j++) {
				prod = table_row[ur1[i+j]];
				lp[j] = prod;
			}
			*lp2 = (*lp2) ^ l;
		}
	}
	free (table_row);
	return;
}

void galois_w08_shift_region_multiply(unsigned char *region,
		int multby,
		int nbytes,
		unsigned char *r2,
		int add)
{
	unsigned char *ur1, *ur2, *ur2top;
	uint16_t product, pp, n1;
	int i, j;
	int num_ones, scratch[8], index, top;

	ur1 = (unsigned char *) region;
	ur2 = (r2 == NULL) ? ur1 : (unsigned char *) r2;
	ur2top = ur2 + nbytes;

	if (multby == 0) {
		if (!add) {
			while (ur2 < ur2top) {
				*ur2 = 0;
				ur2++;
			}
		}
		return;
	}

	pp = prim_poly[8];
	top = 0;
	index = 0;
	num_ones = 0;
	for (i = 0; i < 8; i++) {
		if (multby & (1 << i)) {
			scratch[index] = i;
			index++;
			num_ones++;
			top = i;
		}
	}
	if (!add) {
		for (j = 0; j < nbytes; j++) {
			product = 0;
			n1 = ur1[j];
			for (i = 0; i < num_ones; i++) product ^= (n1 << scratch[i]);
			for (i = top+8; i >= 8; i--) if(product & (1 << i)) product ^= (pp << (i - 8));
			i = (product & 0xff);
			ur2[j] = i;
		}
	} else {
		for (j = 0; j < nbytes; j++) {
			product = 0;
			n1 = ur1[j];
			for (i = 0; i < num_ones; i++) product ^= (n1 << scratch[i]);
			for (i = top+8; i >= 8; i--) if(product & (1 << i)) product ^= (pp << (i - 8));
			i = (product & 0xff);
			ur2[j] ^= i;
		}
	}
	return;
}

void galois_w08_bitmatrix_region_multiply(unsigned char *region,
		int multby,
		int nbytes,
		unsigned char *r2,
		int add)
{
	unsigned char *ur1, *ur2, *ur2top;
	unsigned char prod;
	uint16_t pp;
	int i, j;
	int scratch[8];

	ur1 = (unsigned char *) region;
	ur2 = (r2 == NULL) ? ur1 : (unsigned char *) r2;
	ur2top = ur2 + nbytes;

	if (multby == 0) {
		if (!add) {
			while (ur2 < ur2top) {
				*ur2 = 0;
				ur2++;
			}
		}
		return;
	}

	pp = prim_poly[8];

	for (i = 0; i < 8; i++) {
		scratch[i] = multby;
		multby <<= 1;
		if ((multby & (1 << 8))) {
			multby ^= pp;
		}
	}

	if (!add) {
		for (j = 0; j < nbytes; j++) {
			prod = 0;
			for (i = 0; i < 8; i++) {
				if ((ur1[j] & (1 << i))) prod ^= scratch[i];
			}
			ur2[j] = prod;
		}
	} else {
		for (j = 0; j < nbytes; j++) {
			prod = 0;
			for (i = 0; i < 8; i++) {
				if ((ur1[j] & (1 << i))) prod ^= scratch[i];
			}
			ur2[j] ^= prod;
		}
	}
	return;
	
}

void galois_w08_logtable_region_multiply(unsigned char *region,
		int multby,
		int nbytes,
		unsigned char *r2,
		int add)
{
	unsigned char *ur1, *ur2, *cp;
	int prod;
	int i, log1, j, log2;
	uint64_t l, *lp2, *lptop;
	unsigned char *lp;
	int sol;

	ur1 = (unsigned char *) region;
	ur2 = (r2 == NULL) ? ur1 : (unsigned char *) r2;

	if (multby == 0) {
		if (!add) {
			lp2 = (uint64_t *) ur2;
			ur2 += nbytes;
			lptop = (uint64_t *) ur2;
			while (lp2 < lptop) { *lp2 = 0; lp2++; }
		}
		return;
	}

	if (galois_w08_log_tables[8] == NULL) {
		if (galois_create_log_tables(8) < 0) {
			fprintf(stderr, "galois_w08_logtable_region_multiply -- couldn't make log tables\n");
			exit(1);
		}
	}
	log1 = galois_w08_log_tables[8][multby];

	if (r2 == NULL || !add) {
		for (i = 0; i < nbytes; i++) {
			if (ur1[i] == 0) {
				ur2[i] = 0;
			} else {
				prod = galois_w08_log_tables[8][ur1[i]] + log1;
				ur2[i] = galois_w08_ilog_tables[8][prod];
			}
		}
	} else {
		sol = sizeof(uint64_t);
		lp2 = &l;
		lp = (unsigned char *) lp2;
		for (i = 0; i < nbytes; i += sol) {
			cp = ur2+i;
			lp2 = (uint64_t *) cp;
			for (j = 0; j < sol; j++) {
				if (ur1[i+j] == 0) {
					lp[j] = 0;
				} else {
					log2 = galois_w08_log_tables[8][ur1[i+j]];
					prod = log2 + log1;
					lp[j] = galois_w08_ilog_tables[8][prod];
				}
			}
			*lp2 = (*lp2) ^ l;
		}
	}
	return;
}

void galois_w08_composite_k2_region_multiply (unsigned char *region,
		int multby,
		int nbytes,
		unsigned char *r2,
		int add)
{
	unsigned char *ur1, *ur2, *cp;
	unsigned char prod;
	unsigned char *lp;
	uint64_t l, *lp2, *lptop;
	int sol, i, j;

	uint32_t mask0;
	uint32_t s;
	uint32_t mb0, mb1;
	uint32_t b0, b1;
	uint32_t p2;

	ur1 = region;
	ur2 = (r2 == NULL) ? ur1 : r2;

	if (galois_w08_mult_tables[4] == NULL) {
		if (galois_w08_create_mult_tables(4) < 0) {
			fprintf(stderr, "galois_w08_composite_k2_region_multiply -- couldn't make multiplication tables\n");
			exit(1);
		}
	}

	mask0 = (1 << 4) - 1;
	mb0 = (multby & mask0) << 4;
	mb1 = (multby >> 4) << 4;
	s = 4;

	if (multby == 0) {
		if (!add) {
			lp2 = (uint64_t *) ur2;
			ur2 += nbytes;
			lptop = (uint64_t *) ur2;
			while (lp2 < lptop) { *lp2 = 0; lp2++; }
		}
		return;
	}

	if (r2 == NULL || !add) {
		for (i = 0; i < nbytes; i++) {
			b0 = ur1[i] & mask0;
			b1 = ur1[i] >> 4;
			p2 = galois_w08_mult_tables[4][mb1 | b1];

			prod = galois_w08_mult_tables[4][mb1 | b0];
			prod ^= galois_w08_mult_tables[4][mb0 | b1];
			prod ^= galois_w08_mult_tables[4][(p2 << 4) | s];
			prod <<= 4;
			prod ^= p2;
			prod ^= galois_w08_mult_tables[4][mb0 | b0];

			ur2[i] = prod;
		}
		return;
	} else {
		sol = sizeof(uint64_t);
		lp2 = &l;
		lp = (unsigned char *) lp2;
		for (i = 0; i < nbytes; i += sol) {
			cp = ur2+i;
			lp2 = (uint64_t *) cp;
			for (j = 0; j < sol; j++) {
				b0 = ur1[i+j] & mask0;
				b1 = ur1[i+j] >> 4;
				p2 = galois_w08_mult_tables[4][mb1 | b1];

				prod = galois_w08_mult_tables[4][mb1 | b0];
				prod ^= galois_w08_mult_tables[4][mb0 | b1];
				prod ^= galois_w08_mult_tables[4][(p2 << 4) | s];
				prod <<= 4;
				prod ^= p2;
				prod ^= galois_w08_mult_tables[4][mb0 | b0];

				lp[j] = prod;
			}
			*lp2 = (*lp2) ^ l;
		}
	}
	return;
}

/*  w == 16 implementations:
	log/antilog
	lazy multtable
	shift
	bitmatrix
	composite k == 2
	split composite k == 2 */

void galois_w16_logtable_region_multiply(unsigned char *region,      /* Region to multiply */
		int multby,       /* Number to multiply by */
		int nbytes,        /* Number of bytes in region */
		unsigned char *r2,          /* If r2 != NULL, products go here */
		int add) 
{
	uint16_t *ur1, *ur2, *cp;
	int prod;
	int i, log1, j, log2;
	uint64_t l, *lp2, *lptop;
	uint16_t *lp;
	int sol;

	ur1 = (uint16_t *) region;
	ur2 = (r2 == NULL) ? ur1 : (uint16_t *) r2;
	nbytes /= 2;

	if (multby == 0) {
		if (!add) {
			lp2 = (uint64_t *) ur2;
			ur2 += nbytes;
			lptop = (uint64_t *) ur2;
			while (lp2 < lptop) { *lp2 = 0; lp2++; }
		}
		return;
	}

	if (galois_w16_log_tables[16] == NULL) {
		if (galois_w16_create_log_tables(16) < 0) {
			fprintf(stderr, "galois_w16_logtable_region_multiply -- couldn't make log tables\n");
			exit(1);
		}
	}
	log1 = (int) galois_w16_log_tables[16][multby];

	if (r2 == NULL || !add) {
		for (i = 0; i < nbytes; i++) {
			if (ur1[i] == 0) {
				ur2[i] = 0;
			} else {
				prod = galois_w16_log_tables[16][ur1[i]] + log1;
				ur2[i] = galois_w16_ilog_tables[16][prod];
			}
		}
	} else {
		sol = sizeof(uint64_t)/2;
		lp2 = &l;
		lp = (uint16_t *) lp2;
		for (i = 0; i < nbytes; i += sol) {
			cp = ur2+i;
			lp2 = (uint64_t *) cp;
			for (j = 0; j < sol; j++) {
				if (ur1[i+j] == 0) {
					lp[j] = 0;
				} else {
					log2 = galois_w16_log_tables[16][ur1[i+j]];
					prod = log2 + log1;
					lp[j] = galois_w16_ilog_tables[16][prod];
				}
			}
			*lp2 = (*lp2) ^ l;
		}
	}
	return; 
}

uint16_t *galois_w16_create_lazytable_row (int row_index) {
	uint16_t *row_elts;
	int i;
	int logx;

	if (galois_w16_log_tables[16] == NULL) {
		if (galois_w16_create_log_tables(16) < 0) { return NULL; }
	}

	row_elts = (uint16_t *) malloc(sizeof(uint16_t) * (1 << 16));

	logx = galois_w16_log_tables[16][row_index];
	row_elts[0] = 0;
	for (i = 1; i < (1 << 16); i++) {
		row_elts[i] = galois_w16_ilog_tables[16][logx + galois_w16_log_tables[16][i]];
	}
	return row_elts;
}

void galois_w16_lazytable_region_multiply(unsigned char *region,
		int multby,
		int nbytes,
		unsigned char *r2,
		int add)
{
	uint16_t *table_row;
	uint16_t *ur1, *ur2, *cp;
	uint16_t prod;
	int i, srow, j;
	uint64_t l, *lp2;
	uint16_t *lp;
	int sol;

	ur1 = (uint16_t *) region;
	ur2 = (r2 == NULL) ? ur1 : (uint16_t *) r2;
	nbytes /= 2;

	if (multby == 0) {
		bzero(ur2, nbytes*2);
		return;
	}

	table_row = galois_w16_create_lazytable_row(multby);
	if (table_row == NULL) {
		fprintf(stderr, "galois_16_lazytable_region_multiply -- couldn't make multiplication tables\n");
		exit(1);
	}

	if (r2 == NULL || !add) {
		for (i = 0; i < nbytes; i++) {
			ur2[i] = table_row[ur1[i]];
		}
	} else {
		sol = sizeof(uint64_t)/2;
		lp2 = &l;
		lp = (uint16_t *) lp2;
		for (i = 0; i < nbytes; i += sol) {
			cp = ur2+i;
			lp2 = (uint64_t *) cp;
			for (j = 0; j < sol; j++) {
				lp[j] = table_row[ur1[i+j]];
			}
			*lp2 = (*lp2) ^ l;
		}
	}
	free (table_row);
	return;
}

void galois_w16_shift_region_multiply(unsigned char *region,
		int multby,
		int nbytes,
		unsigned char *r2,
		int add)
{
	uint16_t *ur1, *ur2, *ur2top;
	uint32_t product, pp, n1;
	int i, j;
	int num_ones, scratch[16], index, top;

	ur1 = (uint16_t *) region;
	ur2 = (r2 == NULL) ? ur1 : (uint16_t *) r2;
	nbytes /= sizeof(uint16_t);
	ur2top = ur2 + nbytes;

	if (multby == 0) {
		if (!add) {
			while (ur2 < ur2top) {
				*ur2 = 0;
				ur2++;
			}
		}
		return;
	}

	pp = prim_poly[16];
	top = 0;
	index = 0;
	num_ones = 0;
	for (i = 0; i < 16; i++) {
		if (multby & (1 << i)) {
			scratch[index] = i;
			index++;
			num_ones++;
			top = i;
		}
	}
	if (!add) {
		for (j = 0; j < nbytes; j++) {
			product = 0;
			n1 = ur1[j];
			for (i = 0; i < num_ones; i++) product ^= (n1 << scratch[i]);
			for (i = top+16; i >= 16; i--) if(product & (1 << i)) product ^= (pp << (i - 16));
			i = (product & 0xffff);
			ur2[j] = i;
		}
	} else {
		for (j = 0; j < nbytes; j++) {
			product = 0;
			n1 = ur1[j];
			for (i = 0; i < num_ones; i++) product ^= (n1 << scratch[i]);
			for (i = top+16; i >= 16; i--) if(product & (1 << i)) product ^= (pp << (i - 16));
			i = (product & 0xffff);
			ur2[j] ^= i;
		}
	}
	return;
}

void galois_w16_bitmatrix_region_multiply(unsigned char *region,
		int multby,
		int nbytes,
		unsigned char *r2,
		int add)
{
	uint16_t *ur1, *ur2, *ur2top;
	uint16_t prod;
	int pp;
	int i, j;
	int scratch[16];

	ur1 = (uint16_t *) region;
	ur2 = (r2 == NULL) ? ur1 : (uint16_t *) r2;
	nbytes /= 2;
	ur2top = ur2 + nbytes;

	if (multby == 0) {
		if (!add) {
			while (ur2 < ur2top) {
				*ur2 = 0;
				ur2++;
			}
		}
		return;
	}

	pp = prim_poly[16];

	for (i = 0; i < 16; i++) {
		scratch[i] = multby;
		multby <<= 1;
		if ((multby & (1 << 16))) {
			multby ^= pp;
		}
	}

	if (!add) {
		for (j = 0; j < nbytes; j++) {
			prod = 0;
			for (i = 0; i < 16; i++) {
				if ((ur1[j] & (1 << i))) prod ^= scratch[i];
			}
			ur2[j] = prod;
		}
	} else {
		for (j = 0; j < nbytes; j++) {
			prod = 0;
			for (i = 0; i < 16; i++) {
				if ((ur1[j] & (1 << i))) prod ^= scratch[i];
			}
			ur2[j] ^= prod;
		}
	}
	return;
	
}

void galois_w16_composite_k2_region_multiply (unsigned char *region,
		int multby,
		int nbytes,
		unsigned char *r2,
		int add) 
{
	uint16_t *ur1, *ur2, *cp;
	uint16_t prod;
	uint16_t *lp;
	uint64_t l, *lp2, *lptop;
	int sol, i, j;

	uint16_t mask0;
	uint16_t s;
	uint16_t mb0, mb1;
	uint16_t b0, b1;
	uint16_t p2, u1;


	ur1 = (uint16_t *) region;
	ur2 = (r2 == NULL) ? ur1 : (uint16_t *) r2;
	nbytes /= 2;

	if (multby == 0) {
		if (!add) {
			lp2 = (uint64_t *) ur2;
			ur2 += nbytes;
			lptop = (uint64_t *) ur2;
			while (lp2 < lptop) { *lp2 = 0; lp2++; }
		}
		return;
	}

	if (galois_w08_mult_tables[8] == NULL) {
		if (galois_w08_create_mult_tables(8) < 0) {
			fprintf(stderr, "galois_w16_composite_k2_region_multiply -- couldn't make multiplication tables\n");
			exit(1);
		}
	}

	mask0 = (1 << 8) - 1;
	s = 63;
	mb0 = (multby & mask0) << 8;
	mb1 = (multby >> 8) << 8;

	if (r2 == NULL || !add) {
		for (i = 0; i < nbytes; i++) {
			if (ur1[i] == 0) {
				ur2[i] = 0;
			} else {
				b0 = ur1[i] & mask0;
				b1 = ur1[i] >> 8;
				p2 = galois_w08_mult_tables[8][mb1 | b1];

				prod = galois_w08_mult_tables[8][mb1 | b0];
				prod ^= galois_w08_mult_tables[8][mb0 | b1];
				prod ^= galois_w08_mult_tables[8][(p2 << 8) | s];
				prod <<= 8;
				prod ^= p2;
				prod ^= galois_w08_mult_tables[8][mb0 | b0];

				ur2[i] = prod;
			}
		}
	} else {
		sol = sizeof(uint64_t)/2;
		lp2 = &l;
		lp = (uint16_t *) lp2;
		for (i = 0; i < nbytes; i += sol) {
			cp = &ur2[i];
			lp2 = (uint64_t *) cp;
			for (j = 0; j < sol; j++) {
				u1 = ur1[i+j];
				if (u1 == 0) {
					lp[j] = 0;
				} else {
					b0 = u1 & mask0;
					b1 = u1 >> 8;
					p2 = galois_w08_mult_tables[8][mb1 | b1];

					prod = galois_w08_mult_tables[8][mb1 | b0];
					prod ^= galois_w08_mult_tables[8][mb0 | b1];
					prod ^= galois_w08_mult_tables[8][(p2 << 8) | s];
					prod <<= 8;
					prod ^= p2;
					prod ^= galois_w08_mult_tables[8][mb0 | b0];

					lp[j] = prod;
				}
			}
			*lp2 = (*lp2) ^ l;
		}
	}
	return;
}

void galois_w16_split_composite_k2_region_multiply (unsigned char *region,
		int multby,
		int nbytes,
		unsigned char *r2,
		int add) 
{
	/* low region pointer, high region pointer for destination and source */
	/* multby split into high, low, and s * high */
	int i, j;
	unsigned char *dhp, *shp, *dlp, *slp;
	unsigned char m1, m0, sm1;
	uint32_t mask0;
	unsigned char *tmp;

	/* set region pointers */
	shp = region;
	slp = &region[nbytes/2];
	dhp = (r2 == NULL) ? shp : r2;
	dlp = &dhp[nbytes/2];
	nbytes = nbytes/2;

	/* is multby 0? Speed that up */
	if (multby == 0) {
		if (!add) {
			bzero (dhp, nbytes);
			bzero (dlp, nbytes);
		}
		return;
	}

	/* create w08 mult tables */
	if (galois_w08_mult_tables[8] == NULL) {
		if (galois_w08_create_mult_tables(8) < 0) {
			fprintf(stderr, "galois_w16_split_composite_k2_region_multiply -- couldn't make log tables\n");
			exit(1);
		}
	}
	/* set up constant multipliers */
	mask0 = (1 << 8) -1;
	m1 = (multby >> 8) & mask0;
	m0 = (multby & mask0);
	sm1 = galois_w08_mult_tables[8][(63<<8)|m1];

	if (r2 == NULL) {
		/*****  THIS DOESN'T FUNCTION FOR R2 == NULL AND ADD == 1 *****/
		tmp = (unsigned char *) malloc(sizeof(unsigned char) * nbytes*2);
		dhp = tmp;
		dlp = &tmp[nbytes];
		galois_w08_multtable_region_multiply(slp, m1, nbytes, dhp, add);
		galois_w08_multtable_region_multiply(shp, m0, nbytes, dhp, 1);
		galois_w08_multtable_region_multiply(shp, sm1, nbytes, dhp, 1);

		galois_w08_multtable_region_multiply(shp, m1, nbytes, dlp, add);
		galois_w08_multtable_region_multiply(slp, m0, nbytes, dlp, 1);

		free (region);
		region = tmp;

		return;
	}

	/* w16 region multiplies */
	galois_w08_multtable_region_multiply(slp, m1, nbytes, dhp, add);
	galois_w08_multtable_region_multiply(shp, m0, nbytes, dhp, 1);
	galois_w08_multtable_region_multiply(shp, sm1, nbytes, dhp, 1);

	galois_w08_multtable_region_multiply(shp, m1, nbytes, dlp, add);
	galois_w08_multtable_region_multiply(slp, m0, nbytes, dlp, 1);

	return;
}

/*  w == 32 implementations:
	split w8 tables
	shift
	bitmatrix
	lazy split w8 tables
	lazy split w16 tables
	composite k == 2
	split composite k == 2 */

void galois_w32_split8_table_region_multiply(unsigned char *region,      /* Region to multiply */
		int multby,       /* Number to multiply by */
		int nbytes,        /* Number of bytes in region */
		unsigned char *r2,          /* If r2 != NULL, products go here */
		int add) 
{
	uint32_t *ur1, *ur2, *cp;
	uint64_t *lp2, *lptop;
	int i, j, a, b, i8, j8, k;
	uint32_t acache[4], accumulator;

	ur1 = (uint32_t *) region;
	ur2 = (r2 == NULL) ? ur1 : (uint32_t *) r2;
	nbytes /= sizeof(uint32_t);
	
	if (multby == 0) {
		if (!add) {
			lp2 = (uint64_t *) ur2;
			ur2 += nbytes;
			lptop = (uint64_t *) ur2;
			while (lp2 < lptop) { *lp2 = 0; lp2++; }
		}
		return;
	}

	if (galois_w32_split8_tables[0]== NULL) {
		if (galois_w32_create_split8_tables() < 0) {
			fprintf(stderr, "galois_w32_split8_table_region_multiply -- couldn't make split multiplication tables\n");
			exit(1);
		}
	}

	/* If we're overwriting r2, then we can't do better than just calling split_multiply.
		 We'll inline it here to save on the procedure call overhead */

	i8 = 0;
	for (i = 0; i < 4; i++) {
		acache[i] = (((multby >> i8) & 255) << 8);
		i8 += 8;
	}
	if (!add) {
		for (k = 0; k < nbytes; k++) {
			accumulator = 0;
			for (i = 0; i < 4; i++) {
				a = acache[i];
				j8 = 0;
				for (j = 0; j < 4; j++) {
					b = ((ur1[k] >> j8) & 255);
					accumulator ^= galois_w32_split8_tables[i+j][a|b];
					j8 += 8;
				}
			}
			ur2[k] = accumulator;
		}
	} else {
		for (k = 0; k < nbytes; k++) {
			accumulator = 0;
			for (i = 0; i < 4; i++) {
				a = acache[i];
				j8 = 0;
				for (j = 0; j < 4; j++) {
					b = ((ur1[k] >> j8) & 255);
					accumulator ^= galois_w32_split8_tables[i+j][a|b];
					j8 += 8;
				}
			}
			ur2[k] = (ur2[k] ^ accumulator);
		}
	}
	return;

}

void galois_w32_shift_region_multiply(unsigned char *region,
		int multby,
		int nbytes,
		unsigned char *r2,
		int add)
{
	uint32_t *ur1, *ur2, *ur2top;
	uint64_t product, one, pp, n1;
	int i, j;
	int num_ones, scratch[32], index, top;

	ur1 = (uint32_t *) region;
	ur2 = (r2 == NULL) ? ur1 : (uint32_t *) r2;
	nbytes /= sizeof(uint32_t);
	ur2top = ur2 + nbytes;

	if (multby == 0) {
		if (!add) {
			while (ur2 < ur2top) {
				*ur2 = 0;
				ur2++;
			}
		}
		return;
	}

	one = 1;
	pp = prim_poly[32];
	pp |= (one << 32);
	top = 0;
	index = 0;
	num_ones = 0;
	for (i = 0; i < 32; i++) {
		if (multby & (one << i)) {
			scratch[index] = i;
			index++;
			num_ones++;
			top = i;
		}
	}
	if (!add) {
		for (j = 0; j < nbytes; j++) {
			product = 0;
			n1 = ur1[j];
			for (i = 0; i < num_ones; i++) product ^= (n1 << scratch[i]);
			for (i = top+32; i >= 32; i--) if(product & (one << i)) product ^= (pp << (i - 32));
			i = (product & 0xffffffff);
			ur2[j] = i;
		}
	} else {
		for (j = 0; j < nbytes; j++) {
			product = 0;
			n1 = ur1[j];
			for (i = 0; i < num_ones; i++) product ^= (n1 << scratch[i]);
			for (i = top+32; i >= 32; i--) if(product & (one << i)) product ^= (pp << (i - 32));
			i = (product & 0xffffffff);
			ur2[j] ^= i;
		}
	}
	return;
}

void galois_w32_bitmatrix_region_multiply(unsigned char *region,
		int multby,
		int nbytes,
		unsigned char *r2,
		int add)
{
	uint32_t *ur1, *ur2, *ur2top;
	uint32_t prod;
	uint64_t pp;
	int i, j;
	int scratch[32];

	ur1 = (uint32_t *) region;
	ur2 = (r2 == NULL) ? ur1 : (uint32_t *) r2;
	nbytes /= 4;
	ur2top = ur2 + nbytes;

	if (multby == 0) {
		if (!add) {
			while (ur2 < ur2top) {
				*ur2 = 0;
				ur2++;
			}
		}
		return;
	}

	pp = prim_poly[32];

	for (i = 0; i < 32; i++) {
		scratch[i] = multby;
		if (multby & (1 << 31)) {
			multby <<= 1;
			multby ^= pp;
		} else {
			multby <<= 1;
		}
	}

	if (!add) {
		for (j = 0; j < nbytes; j++) {
			prod = 0;
			for (i = 0; i < 32; i++) {
				if ((ur1[j] & (1 << i))) prod ^= scratch[i];
			}
			ur2[j] = prod;
		}
	} else {
		for (j = 0; j < nbytes; j++) {
			prod = 0;
			for (i = 0; i < 32; i++) {
				if ((ur1[j] & (1 << i))) prod ^= scratch[i];
			}
			ur2[j] ^= prod;
		}
	}
	return;
	
}

uint32_t **galois_w32_create_lazy_split8_row(int multby) {
	int p1, p2, index, ishift, jshift, i, j, i8, j8, k;
	uint32_t p1elt, p2elt, acache[4], **table_row;

	table_row = (uint32_t **) malloc(sizeof(uint32_t *) * 16);
	if (table_row == NULL) { return NULL; }

	for (i = 0; i < 16; i++) {
		table_row[i] = (uint32_t *) malloc(sizeof(uint32_t) * (1 << 8));
		if (table_row[i] == NULL) {
			for (i--; i >= 0; i--) free(table_row[i]);
			free(table_row);
			return NULL;
		}
	}

	i8 = 0;
	for (i = 0; i < 4; i++) {
		acache[i] = (multby >> i8) & 255;
		i8 += 8;
	}

	for (i = 0; i < 4; i++) {
		for (j = 0; j < 4; j++) {
			if (i+j > 3) {
				jshift = (i+j-3) * 8;
				ishift = 3 * 8;
			} else {
				jshift = (i+j) * 8;
				ishift = 0;
			}

			index = 0;
			p1elt = (acache[i] << ishift);
			for (p2 = 0; p2 < 256; p2++) {
				p2elt = (p2 << jshift);
				table_row[i*4+j][index] = (uint32_t) galois_bitmatrix_multiply(p1elt, p2elt, 32);
				index++;
			}
		}
	}
	return table_row;
}

void galois_w32_lazy_split8_table_region_multiply(unsigned char *region,
		int multby,
		int nbytes,
		unsigned char *r2,
		int add)
{
	uint32_t *ur1, *ur2;
	int i, j, b, i8, j8, k;
	int p1, p2;
	uint32_t accumulator, *table;

	uint32_t **table_row;

	ur1 = (uint32_t *) region;
	ur2 = (r2 == NULL) ? ur1 : (uint32_t *) r2;
	nbytes /= sizeof(uint32_t);

	if (multby == 0) {
		if (!add) {
			bzero(r2, nbytes * sizeof(uint32_t));
		}
		return;
	}

	table_row = galois_w32_create_lazy_split8_row(multby);

	if (!add) {
		for (k = 0; k < nbytes; k++) {
			accumulator = 0;
			for (i = 0; i < 4; i++) {
				j8 = 0;
				for (j = 0; j < 4; j++) {
					b = ((ur1[k] >> j8) & 255);
					accumulator ^= table_row[i*4 + j][b];
					j8 += 8;
				}
			}
			ur2[k] = accumulator;
		}
	} else {
		for (k = 0; k < nbytes; k++) {
			accumulator = 0;
			j8 = 0;
			for (j = 0; j < 4; j++) {
				b = (ur1[k] >> j8) & 255;
				for (i = 0; i < 4; i++) {
					accumulator ^= table_row[i*4 + j][b];
				}
				j8 += 8;
			}
			ur2[k] = (ur2[k] ^ accumulator);
		}
	}
	for (i = 0; i < 16; i++) free(table_row[i]);
	free(table_row);
	return;
}

uint32_t **galois_w32_create_lazy_split16_row(int multby) {
	int p1, p2, index, ishift, jshift, i, j, i16, k;
	uint32_t p1elt, p2elt, acache[2], **table_row;

	table_row = (uint32_t **) malloc(sizeof(uint32_t *) * 4);
	if (table_row == NULL) { return NULL; }

	for (i = 0; i < 4; i++) {
		table_row[i] = (uint32_t *) malloc(sizeof(uint32_t) * (1 << 16));
		if (table_row[i] == NULL) {
			for (i--; i >= 0; i--) free(table_row[i]);
			free(table_row);
			return NULL;
		}
	}

	i16 = 0;
	for (i = 0; i < 2; i++) {
		acache[i] = (multby >> i16) & 65535;
		i16 += 16;
	}

	for (i = 0; i < 2; i++) {
		for (j = 0; j < 2; j++) {
			if (i+j > 1) {
				jshift = (i+j-1) * 16;
				ishift = 1 * 16;
			} else {
				jshift = (i+j) * 16;
				ishift = 0;
			}

			index = 0;
			p1elt = (acache[i] << ishift);
			for (p2 = 0; p2 < 65536; p2++) {
				p2elt = (p2 << jshift);
				table_row[i*2+j][index] = (uint32_t) galois_bitmatrix_multiply(p1elt, p2elt, 32);
				index++;
			}
		}
	}
	return table_row;
}

void galois_w32_lazy_split16_table_region_multiply(unsigned char *region,
		int multby,
		int nbytes,
		unsigned char *r2,
		int add)
{
	uint32_t *ur1, *ur2;
	int i, j, b, j16, k;
	int p1, p2;
	uint32_t accumulator, *table;

	uint32_t **table_row;

	ur1 = (uint32_t *) region;
	ur2 = (r2 == NULL) ? ur1 : (uint32_t *) r2;
	nbytes /= sizeof(uint32_t);

	if (multby == 0) {
		if (!add) {
			bzero(r2, nbytes * sizeof(uint32_t));
		}
		return;
	}

	table_row = galois_w32_create_lazy_split16_row(multby);

	if (!add) {
		for (k = 0; k < nbytes; k++) {
			accumulator = 0;
			for (i = 0; i < 2; i++) {
				j16 = 0;
				for (j = 0; j < 2; j++) {
					b = ((ur1[k] >> j16) & 65535);
					accumulator ^= table_row[i*2 + j][b];
					j16 += 16;
				}
			}
			ur2[k] = accumulator;
		}
	} else {
		for (k = 0; k < nbytes; k++) {
			accumulator = 0;
			j16 = 0;
			for (j = 0; j < 2; j++) {
				b = (ur1[k] >> j16) & 65535;
				for (i = 0; i < 2; i++) {
					accumulator ^= table_row[i*2 + j][b];
				}
				j16 += 16;
			}
			ur2[k] = (ur2[k] ^ accumulator);
		}
	}
	for (i = 0; i < 4; i++) free(table_row[i]);
	free(table_row);
	return;
	
}

/* This one needs some work, as well.  Try to inline.  Pay attention to unsigned integers not acting as indices for ilog tables?  But I'm only using mult.... */
void galois_w32_composite_k2_region_multiply (unsigned char *region,
		int multby,
		int nbytes,
		unsigned char *r2,
		int add) 
{
	uint32_t *ur1, *ur2, *cp, *lp, l, *lp2;
	uint32_t i, j;

	int mblog1, mblog0, blog1, blog0, sum_log;
	uint32_t mb;
	uint32_t prod;
	uint32_t mask0, s, p2, p3;

	mb = multby;
	ur1 = (uint32_t *) region;
	ur2 = (r2 == NULL) ? ur1 : (uint32_t *) r2;
	nbytes /= sizeof(uint32_t);

	if (multby == 0) {
		if (!add) {
			lp = (uint32_t *) ur2;
			ur2 += nbytes;
			lp2 = (uint32_t *) ur2;
			while (lp < lp2) { *lp = 0; lp++; }
		}
		return;
	}

	if (galois_w16_log_tables[16] == NULL) {
		if (galois_w16_create_log_tables(16) < 0) {
			fprintf(stderr, "galois_w32_composite_k2_region_multiply -- couldn't make log tables\n");
			exit(1);
		}
	}

	s = 40188;
	mask0 = (1 << 16) - 1;
	mblog1 = galois_w16_log_tables[16][((mb >> 16) & mask0)];
	mblog0 = galois_w16_log_tables[16][(mb & mask0)];

	/* This has not been inlined because logtable multiplication requires both factors zero checked.
	   Direct inlining would be messy (5 if-blocks) and the other alternative requires 8 if blocks
	   to check the separate values in multby and ur1[i] */
	if (r2 == NULL || !add) {
		for (i = 0; i < nbytes; i++) {
			if ((int) ur1[i] == 0) {
				ur2[i] = 0;
			} else {
								
				p2 = galois_w16_logtable_multiply(mb >> 16, ur1[i] >> 16, 16);
				p3 = p2;
				p2 = galois_w16_logtable_multiply(p2, s, 16);

				prod = galois_w16_logtable_multiply(mb >> 16, ur1[i] & mask0, 16);
				prod ^= galois_w16_logtable_multiply(mb & mask0, ur1[i] >> 16, 16);
				prod ^= p2;
				prod <<= 16;
				prod ^= p3;
				prod ^= galois_w16_logtable_multiply(mb & mask0, ur1[i] & mask0, 16);

				ur2[i] = prod;
			}
		}
	} else {
		lp = (uint32_t *) &l;
		for (i = 0; i < nbytes; i ++) {
			cp = ur2+i;
			lp2 = (uint32_t *) cp;
			if (ur1[i] == 0) {
				l = 0;
			} else {
				p2 = galois_w16_logtable_multiply(mb >> 16, ur1[i] >> 16, 16);

				prod = galois_w16_logtable_multiply(mb >> 16, ur1[i] & mask0, 16);
				prod ^= galois_w16_logtable_multiply(mb & mask0, ur1[i] >> 16, 16);
				prod ^= galois_w16_logtable_multiply(s, p2, 16);
				prod <<= 16;
				prod ^= p2;
				prod ^= galois_w16_logtable_multiply(mb & mask0, ur1[i] & mask0, 16);

				l = prod;
			}
			*lp2 = (*lp2) ^ l;
		}
	}
	return;
}

void galois_w32_split_composite_k2_region_multiply (unsigned char *region,
		int multby,
		int nbytes,
		unsigned char *r2,
		int add) 
{
	/* low region pointer, high region pointer for destination and source */
	/* multby split into high, low, and s * high */
	int i, j;
	unsigned char *dhp, *shp, *dlp, *slp;
	uint16_t m1, m0, sm1;
	uint32_t mask0;
	unsigned char *tmp;

	/* set region pointers */
	shp = region;
	slp = &region[nbytes/2];
	dhp = (r2 == NULL) ? shp : r2;
	dlp = &dhp[nbytes/2];
	nbytes = nbytes/2;

	/* is multby 0? Speed that up */
	if (multby == 0) {
		if (!add) {
			bzero (dhp, nbytes);
			bzero (dlp, nbytes);
		}
		return;
	}

	/* create w16 log tables */
	if (galois_w16_log_tables[16] == NULL) {
		if (galois_w16_create_log_tables(16) < 0) {
			fprintf(stderr, "galois_w32_split_composite_k2_region_multiply -- couldn't make log tables\n");
			exit(1);
		}
	}

	/* set up constant multipliers */
	mask0 = (1 << 16) -1;
	m1 = (multby >> 16) & mask0;
	m0 = (multby & mask0);

	sm1 = galois_w16_logtable_multiply(m1, 40188, 16);

	if (r2 == NULL) {
		tmp = (unsigned char *) malloc(sizeof(unsigned char) * nbytes*2);
		dhp = tmp;
		dlp = &tmp[nbytes];
		galois_w16_logtable_region_multiply(slp, m1, nbytes, dhp, add);
		galois_w16_logtable_region_multiply(shp, m0, nbytes, dhp, 1);
		galois_w16_logtable_region_multiply(shp, sm1, nbytes, dhp, 1);

		galois_w16_logtable_region_multiply(shp, m1, nbytes, dlp, add);
		galois_w16_logtable_region_multiply(slp, m0, nbytes, dlp, 1);

		free (region);
		region = tmp;

		return;
	}

	/* w16 region multiplies */
	galois_w16_logtable_region_multiply(slp, m1, nbytes, dhp, add);
	galois_w16_logtable_region_multiply(shp, m0, nbytes, dhp, 1);
	galois_w16_logtable_region_multiply(shp, sm1, nbytes, dhp, 1);

	galois_w16_logtable_region_multiply(shp, m1, nbytes, dlp, add);
	galois_w16_logtable_region_multiply(slp, m0, nbytes, dlp, 1);

	return;
}

/*  w == 64 implementations:
	split w8 tables
	lazy split w8 tables
	lazy split w16 tables
	shift
	bitmatrix
	composite k == 2  (using split w8 for l == 32)
	split composite k == 2 
	recursive composite k == 2 (calls split composite k==2) */

void galois_w64_split8_table_region_multiply (unsigned char *region,
		uint64_t multby,
		int nbytes,
		unsigned char *r2,
		int add)
{
	uint64_t *ur1, *ur2, *cp;
	uint64_t *lp2, *lptop;
	uint64_t i, j, a, b, accumulator, i8, j8, k;
	uint64_t acache[8];

	ur1 = (uint64_t *) region;
	ur2 = (r2 == NULL) ? ur1 : (uint64_t *) r2;
	nbytes /= sizeof(uint64_t);

	if (multby == 0) {
		if (!add) {
			lp2 = (uint64_t *) ur2;
			ur2 += nbytes;
			lptop = (uint64_t *) ur2;
			while (lp2 < lptop) { *lp2 = 0; lp2++; }
		}
		return;
	}

	if (galois_w64_split8_tables[0]== NULL) {
		if (galois_w64_create_split8_tables() < 0) {
			fprintf(stderr, "galois_64_region_multiply -- couldn't make split multiplication tables\n");
			exit(1);
		}
	}

	/* If we're overwriting r2, then we can't do better than just calling split_multiply.
		 We'll inline it here to save on the procedure call overhead */

	i8 = 0;
	for (i = 0; i < 8; i++) {
		acache[i] = (((multby >> i8) & 255) << 8);
		i8 += 8;
	}
	if (!add) {
		for (k = 0; k < nbytes; k++) {
			accumulator = 0;
			for (i = 0; i < 8; i++) {
				a = acache[i];
				j8 = 0;
				for (j = 0; j < 8; j++) {
					b = ((ur1[k] >> j8) & 255);
					accumulator ^= galois_w64_split8_tables[i+j][a|b];
					j8 += 8;
				}
			}
			ur2[k] = accumulator;
		}
	} else {
		for (k = 0; k < nbytes; k++) {
			accumulator = 0;
			for (i = 0; i < 8; i++) {
				a = acache[i];
				j8 = 0;
				for (j = 0; j < 8; j++) {
					b = ((ur1[k] >> j8) & 255);
					accumulator ^= galois_w64_split8_tables[i+j][a|b];
					j8 += 8;
				}
			}
			ur2[k] = (ur2[k] ^ accumulator);
		}
	}
	return;
}

uint64_t **galois_w64_create_lazy_split8_row(uint64_t multby) {
	uint64_t p1, p2, index, ishift, jshift, i, j, i8, j8, k;
	uint64_t p1elt, p2elt, acache[8], **table_row;

	table_row = (uint64_t **) malloc(sizeof(uint64_t *) * 64);
	if (table_row == NULL) { return NULL; }

	for (i = 0; i < 64; i++) {
		table_row[i] = (uint64_t *) malloc(sizeof(uint64_t) * (1 << 8));
		if (table_row[i] == NULL) {
			for (i--; i >= 0; i--) free(table_row[i]);
			free(table_row);
			return NULL;
		}
	}

	i8 = 0;
	for (i = 0; i < 8; i++) {
		acache[i] = (multby >> i8) & 255;
		i8 += 8;
	}

	for (i = 0; i < 8; i++) {
		for (j = 0; j < 8; j++) {
			if (i+j > 7) {
				jshift = (i+j-7) * 8;
				ishift = 7 * 8;
			} else {
				jshift = (i+j) * 8;
				ishift = 0;
			}

			index = 0;
			p1elt = (acache[i] << ishift);
			for (p2 = 0; p2 < 256; p2++) {
				p2elt = (p2 << jshift);
				table_row[i*8+j][index] = (uint64_t) galois_w64_bitmatrix_multiply(p1elt, p2elt);
				index++;
			}
		}
	}
	return table_row;
}

void galois_w64_lazy_split8_table_region_multiply (unsigned char *region,
		uint64_t multby,
		int nbytes,
		unsigned char *r2,
		int add)
{
	uint64_t *ur1, *ur2, *cp, *ur2top;
	uint64_t *lp2, *lptop;
	uint64_t i, j, b, accumulator, i8, j8, k;
	uint64_t **table_row;

	ur1 = (uint64_t *) region;
	ur2 = (r2 == NULL) ? ur1 : (uint64_t *) r2;
	nbytes /= sizeof(uint64_t);
	ur2top = ur2 + nbytes;

	if (multby == 0) {
		if (!add) {
			bzero(r2, nbytes * sizeof(uint64_t));
		}
		return;
	}

	/* If we're overwriting r2, then we can't do better than just calling split_multiply.
		 We'll inline it here to save on the procedure call overhead */

	table_row = galois_w64_create_lazy_split8_row(multby);
	i8 = 0;
	if (!add) {
		for (k = 0; k < nbytes; k++) {
			accumulator = 0;
			for (i = 0; i < 8; i++) {
				j8 = 0;
				for (j = 0; j < 8; j++) {
					b = ((ur1[k] >> j8) & 255);
					accumulator ^= table_row[i*8 + j][b];
					j8 += 8;
				}
			}
			ur2[k] = accumulator;
		}
	} else {
		for (k = 0; k < nbytes; k++) {
			accumulator = 0;
			for (i = 0; i < 8; i++) {
				j8 = 0;
				for (j = 0; j < 8; j++) {
					b = ((ur1[k] >> j8) & 255);
					accumulator ^= table_row[i*8 + j][b];
					j8 += 8;
				}
			}
			ur2[k] = (ur2[k] ^ accumulator);
		}
	}
	for (i = 0; i < 16; i++) free(table_row[i]);
	free(table_row);
	return;
}

uint64_t **galois_w64_create_lazy_split16_row(uint64_t multby) {
	int index, ishift, jshift, i, j, i16;
	uint64_t p2, p1elt, p2elt, acache[4], **table_row;

	table_row = (uint64_t **) malloc(sizeof(uint64_t *) * 16);
	if (table_row == NULL) { return NULL; }

	for (i = 0; i < 16; i++) {
		table_row[i] = (uint64_t *) malloc(sizeof(uint64_t) * (1 << 16));
		if (table_row[i] == NULL) {
			for (i--; i >= 0; i--) free(table_row[i]);
			free(table_row);
			return NULL;
		}
	}

	i16 = 0;
	for (i = 0; i < 4; i++) {
		acache[i] = (rshift(multby, i16) & 65535);
		i16 += 16;
	}

	for (i = 0; i < 4; i++) {
		for (j = 0; j < 4; j++) {
			if (i+j > 3) {
				jshift = (i+j-3) * 16;
				ishift = 3 * 16;
			} else {
				jshift = (i+j) * 16;
				ishift = 0;
			}

			index = 0;
			p1elt = lshift(acache[i], ishift);
			for (p2 = 0; p2 < 65536; p2++) {
				p2elt = lshift(p2, jshift);
				table_row[i*4+j][index] = galois_w64_bitmatrix_multiply(p1elt, p2elt);
				index++;
			}
		}
	}
	return table_row;
}

void galois_w64_lazy_split16_table_region_multiply (unsigned char *region,
		uint64_t multby,
		int nbytes,
		unsigned char *r2,
		int add)
{
	uint64_t *ur1, *ur2;
	int i, j, b, j16, k;
	int p1, p2;
	uint64_t accumulator;

	uint64_t **table_row;

	ur1 = (uint64_t *) region;
	ur2 = (r2 == NULL) ? ur1 : (uint64_t *) r2;
	nbytes /= sizeof(uint64_t);

	if (multby == 0) {
		if (!add) {
			bzero(r2, nbytes * sizeof(uint64_t));
		}
		return;
	}

	table_row = galois_w64_create_lazy_split16_row(multby);

	if (!add) {
		for (k = 0; k < nbytes; k++) {
			accumulator = 0;
			for (i = 0; i < 4; i++) {
				j16 = 0;
				for (j = 0; j < 4; j++) {
					b = (rshift(ur1[k], j16) & 65535);
					accumulator ^= table_row[i*4 + j][b];
					j16 += 16;
				}
			}
			ur2[k] = accumulator;
		}
	} else {
		for (k = 0; k < nbytes; k++) {
			accumulator = 0;
			for (i = 0; i < 4; i++) {
				j16 = 0;
				for (j = 0; j < 4; j++) {
					b = (rshift(ur1[k], j16) & 65535);
					accumulator ^= table_row[i*4 + j][b];
					j16 += 16;
				}
			}
			ur2[k] ^= accumulator;
		}
	}
	for (i = 0; i < 16; i++) free(table_row[i]);
	free(table_row);
	return;
	
}

void galois_w64_shift_region_multiply (unsigned char *region,
	uint64_t multby,
	int nbytes,
	unsigned char *r2,
	int add)
{
	uint64_t *ur1, *ur2, *ur2top;
	uint64_t product[2], one, pp;
	int i, j;
	int num_ones, scratch[64], index, top;

	ur1 = (uint64_t *) region;
	ur2 = (r2 == NULL) ? ur1 : (uint64_t *) r2;
	nbytes /= sizeof(uint64_t);
	ur2top = ur2 + nbytes;

	if (multby == 0) {
		if (!add) {
			while (ur2 < ur2top) {
				*ur2 = 0;
				ur2++;
			}
		}
		return;
	}

	one = 1;
	pp = 0x1b;
	top = 0;
	index = 0;
	num_ones = 0;
	for (i = 0; i < 64; i++) {
		if (multby & (one << i)) {
			scratch[index] = i;
			index++;
			num_ones++;
			top = i;
		}
	}
	if (!add) {
		for (j = 0; j < nbytes; j++) {
			product[0] = 0;
			product[1] = 0;
			for (i = 0; i < num_ones; i++) {
				product[0] ^= rshift(ur1[j], (64-scratch[i]));
				product[1] ^= lshift(ur1[j], scratch[i]);
				//product[0] ^= (ur1[j] >> (64 - scratch[i]));
				//product[1] ^= (ur1[j] << scratch[i]);
			}
			for (i = top+64; i >= 64; i--) {
				if(product[0] & lshift(one, i-64)) {
					product[0] ^= rshift(pp, 128-i);
					product[1] ^= lshift(pp, i-64);
					//product[0] ^= (pp >> (128 - i));
					//product[1] ^= (pp << (i-64));
				}
			}
			ur2[j] = product[1];
		}
	} else {
		for (j = 0; j < nbytes; j++) {
			product[0] = 0;
			product[1] = 0;
			for (i = 0; i < num_ones; i++) {
				product[0] ^= rshift(ur1[j], (64-scratch[i]));
				product[1] ^= lshift(ur1[j], scratch[i]);
				//product[0] ^= (ur1[j] >> (64 - scratch[i]));
				//product[1] ^= (ur1[j] << scratch[i]);
			}
			for (i = top+64; i >= 64; i--) {
				if(product[0] & lshift(one, i-64)) {
					product[0] ^= rshift(pp, 128-i);
					product[1] ^= lshift(pp, i-64);
					//product[0] ^= (pp >> (128 - i));
					//product[1] ^= (pp << (i-64));
				}
			}
			ur2[j] ^= product[1];
		}
	}
	return;

}

void galois_w64_bitmatrix_region_multiply (unsigned char *region,
		uint64_t multby,
		int nbytes,
		unsigned char *r2,
		int add)
{
	uint64_t *ur1, *ur2, *ur2top;
	uint64_t pp, one, prod;
	int i, j;
	int scratch[64];

	ur1 = (uint64_t *) region;
	ur2 = (r2 == NULL) ? ur1 : (uint64_t *) r2;
	nbytes /= 8;
	ur2top = ur2 + nbytes;

	if (multby == 0) {
		if (!add) {
			while (ur2 < ur2top) {
				*ur2 = 0;
				ur2++;
			}
		}
		return;
	}

	one = 1;
	pp = 0x1b;

	for (i = 0; i < 64; i++) {
		scratch[i] = multby;
		if (multby & (one << 63)) {
			multby <<= 1;
			multby ^= pp;
		} else {
			multby <<= 1;
		}
	}

	if (!add) {
		for (j = 0; j < nbytes; j++) {
			prod = 0;
			for (i = 0; i < 64; i++) {
				if ((ur1[j] & (one << i))) prod ^= scratch[i];
			}
			ur2[j] = prod;
		}
	} else {
		for (j = 0; j < nbytes; j++) {
			prod = 0;
			for (i = 0; i < 64; i++) {
				if ((ur1[j] & (one << i))) prod ^= scratch[i];
			}
			ur2[j] ^= prod;
		}
	}
	return;

}

void galois_w64_composite_k2_region_multiply (unsigned char *region,
		uint64_t multby,
		int nbytes,
		unsigned char *r2,
		int add)
{
	uint64_t *ur1, *ur2, *lp, *lp2;
	int i, l;

	uint64_t mb;
	uint64_t prod;
	uint64_t mask0, s, p2, p3;
	uint64_t one;

	mb = multby;
	ur1 = (uint64_t *) region;
	ur2 = (r2 == NULL) ? ur1 : (uint64_t *) r2;
	nbytes /= 8;

	if (multby == 0) {
		if (!add) {
			lp = (uint64_t *) ur2;
			ur2 += nbytes;
			lp2 = (uint64_t *) ur2;
			while (lp < lp2) { *lp = 0; lp++; }
		}
		return;
	}

	if (galois_w32_split8_tables[0] == NULL) {
		if (galois_w32_create_split8_tables() < 0) {
			fprintf(stderr, "galois_w64_composite_k2_region_multiply: Error creating tables for w=32\n");
			exit(1);
		}
	}

	s = 2;
	l = 32;
	one = 1;
	mask0 = lshift(one, l) - 1;

	if (r2 == NULL || !add) {
		for (i = 0; i < nbytes; i++) {

			p2 = galois_w32_split8_table_multiply(((multby >> l) & mask0), ((ur1[i] >> l) & mask0));
			p3 = p2;
			p2 = galois_w32_split8_table_multiply((s & mask0), (p2 & mask0));

			prod = galois_w32_split8_table_multiply(((multby >> l) & mask0), (ur1[i] & mask0));
			prod ^= galois_w32_split8_table_multiply((multby & mask0), ((ur1[i] >> l) & mask0));
			prod ^= p2;
			prod <<= l;
			prod ^= p3;
			prod ^= galois_w32_split8_table_multiply((multby & mask0), (ur1[i] & mask0));
			
			ur2[i] = prod;
		}
	} else {
		for (i = 0; i < nbytes; i++) {
			lp2 = (uint64_t *) ur2 + i;
			
			p2 = galois_w32_split8_table_multiply(((multby >> l) & mask0), ((ur1[i] >> l) & mask0));
			p3 = p2;
			p2 = galois_w32_split8_table_multiply((s & mask0), (p2 & mask0));

			prod = galois_w32_split8_table_multiply(((multby >> l) & mask0), (ur1[i] & mask0));
			prod ^= galois_w32_split8_table_multiply((multby & mask0), ((ur1[i] >> l) & mask0));
			prod ^= p2;
			prod <<= l;
			prod ^= p3;
			prod ^= galois_w32_split8_table_multiply((multby & mask0), (ur1[i] & mask0));

			*lp2 = (*lp2) ^ prod;
		}
	}
	return;
	
}
		
void galois_w64_split_composite_k2_region_multiply (unsigned char *region,
		uint64_t multby,
		int nbytes,
		unsigned char *r2,
		int add)
{
	/* low region pointer, high region pointer for destination and source */
	/* multby split into high, low, and s * high */
	int i, j;
	unsigned char *dhp, *shp, *dlp, *slp;
	uint64_t m1, m0, sm1;
	uint64_t mask0, one;
	unsigned char *tmp;

	/* set region pointers */
	shp = region;
	slp = region + nbytes/2;
	dhp = (r2 == NULL) ? shp : r2;
	dlp = dhp + nbytes/2;
	nbytes = nbytes/2;

	/* is multby 0? Speed that up */
	if (multby == 0) {
		if (!add) {
			bzero (dhp, nbytes);
			bzero (dlp, nbytes);
		}
		return;
	}

	/* create w32 split tables */
	if (galois_w32_split8_tables[0] == NULL) {
		if (galois_w32_create_split8_tables() < 0) {
			fprintf(stderr, "galois_w64_split_composite_k2_region_multiply -- couldn't make split tables\n");
			exit(1);
		}
	}

	/* set up constant multipliers */
	one = 1;
	mask0 = 0xffffffff;
	m1 = rshift(multby, 32) & mask0;
	m0 = (multby & mask0);

	sm1 = galois_w32_split8_table_multiply(m1, 2);

	if (r2 == NULL) {
		tmp = (unsigned char *) malloc(sizeof(unsigned char) * nbytes*2);
		dhp = tmp;
		dlp = tmp + nbytes;
		galois_w32_split8_table_region_multiply(slp, m1, nbytes, dhp, add);
		galois_w32_split8_table_region_multiply(shp, m0, nbytes, dhp, 1);
		galois_w32_split8_table_region_multiply(shp, sm1, nbytes, dhp, 1);

		galois_w32_split8_table_region_multiply(shp, m1, nbytes, dlp, add);
		galois_w32_split8_table_region_multiply(slp, m0, nbytes, dlp, 1);

		free (region);
		region = tmp;

		return;
	}

	/* w32 region multiplies */
	galois_w32_split8_table_region_multiply(slp, m1, nbytes, dhp, add);
	galois_w32_split8_table_region_multiply(shp, m0, nbytes, dhp, 1);
	galois_w32_split8_table_region_multiply(shp, sm1, nbytes, dhp, 1);

	galois_w32_split8_table_region_multiply(shp, m1, nbytes, dlp, add);
	galois_w32_split8_table_region_multiply(slp, m0, nbytes, dlp, 1);

	return;

}

void galois_w64_recursive_composite_k2_region_multiply (unsigned char *region,
		uint64_t multby,
		int nbytes,
		unsigned char *r2,
		int add)
{
	/* low region pointer, high region pointer for destination and source */
	/* multby split into high, low, and s * high */
	int i, j;
	unsigned char *dhp, *shp, *dlp, *slp;
	uint64_t m1, m0, sm1;
	uint64_t mask0, one;
	unsigned char *tmp;

	/* set region pointers */
	shp = region;
	slp = region + nbytes/2;
	dhp = (r2 == NULL) ? shp : r2;
	dlp = dhp + nbytes/2;
	nbytes = nbytes/2;

	/* is multby 0? Speed that up */
	if (multby == 0) {
		if (!add) {
			bzero (dhp, nbytes);
			bzero (dlp, nbytes);
		}
		return;
	}

	/* set up constant multipliers */
	one = 1;
	mask0 = 0xffffffff;
	m1 = rshift(multby, 32) & mask0;
	m0 = (multby & mask0);

	sm1 = galois_w32_composite_k2_multiply(m1, 2);

	if (r2 == NULL) {
		tmp = (unsigned char *) malloc(sizeof(unsigned char) * nbytes*2);
		dhp = tmp;
		dlp = tmp + nbytes;
		galois_w32_split_composite_k2_region_multiply(slp, m1, nbytes, dhp, add);
		galois_w32_split_composite_k2_region_multiply(shp, m0, nbytes, dhp, 1);
		galois_w32_split_composite_k2_region_multiply(shp, sm1, nbytes, dhp, 1);

		galois_w32_split_composite_k2_region_multiply(shp, m1, nbytes, dlp, add);
		galois_w32_split_composite_k2_region_multiply(slp, m0, nbytes, dlp, 1);

		free (region);
		region = tmp;

		return;
	}

	/* w32 region multiplies */
	galois_w32_split_composite_k2_region_multiply(slp, m1, nbytes, dhp, add);
	galois_w32_split_composite_k2_region_multiply(shp, m0, nbytes, dhp, 1);
	galois_w32_split_composite_k2_region_multiply(shp, sm1, nbytes, dhp, 1);

	galois_w32_split_composite_k2_region_multiply(shp, m1, nbytes, dlp, add);
	galois_w32_split_composite_k2_region_multiply(slp, m0, nbytes, dlp, 1);

	return;

}

/*  w == 128
	shift
	bitmatrix
	split w8 tables
	lazy split w8 tables
	composite k == 2 (using split table for l == 64)
	split composite k == 2
	recursive composite k == 2 (calling galois_w64_recursive_composite_k2_RM) */

void galois_w128_shift_region_multiply (unsigned char *region,
		uint64_t multby[2],
		int nbytes,
		unsigned char *r2,
		int add) 
{
	uint64_t *ur1, *ur2, *ur2top;
	uint64_t product[4], one, pp;
	int i, j;
	int num_ones, scratch[128], index, top;

	ur1 = (uint64_t *) region;
	ur2 = (r2 == NULL) ? ur1 : (uint64_t *) r2;
	nbytes /= sizeof(uint64_t);
	ur2top = ur2 + nbytes;

	if (multby == 0) {
		if (!add) {
			while (ur2 < ur2top) {
				*ur2 = 0;
				ur2++;
			}
		}
		return;
	}

	one = 1;
	pp = 0x87;
	top = 0;
	index = 0;
	num_ones = 0;
	for (i = 0; i < 64; i++) {
		if ((multby[1] & lshift(one, i))) {
			scratch[index] = i;
			index++;
			num_ones++;
			if (i > top) top = i;
		}
		if ((multby[0] & lshift(one, i))) {
			scratch[index] = i+64;
			index++;
			num_ones++;
			top = i+64;
		}
	}
	if (!add) {
		for (j = 0; j < nbytes; j+=2) {
			product[0] = 0;
			product[1] = 0;
			product[2] = 0;
			product[3] = 0;
			for (i = 0; i < num_ones; i++) {
				if (scratch[i]<64) {
					product[1] ^= rshift(ur1[j], 64-scratch[i]);
					product[2] ^= lshift(ur1[j], scratch[i]);
					product[2] ^= rshift(ur1[j+1], 64-scratch[i]);
					product[3] ^= lshift(ur1[j+1], scratch[i]);
				} else if (scratch[i] >= 64) {
					product[0] ^= rshift(ur1[j], 128-scratch[i]);
					product[1] ^= lshift(ur1[j], scratch[i]-64);
					product[1] ^= rshift(ur1[j+1], 128-scratch[i]);
					product[2] ^= lshift(ur1[j+1], scratch[i]-64);
				}
			}
			for (i = top+128; i >= 128; i--) {
				if ((i >= 192) && (product[0] & (lshift(one, i-192)))) {
					product[1] ^= rshift(pp, 64-(i-192));
					product[2] ^= lshift(pp, i-192);
				} else if ((i >= 128) && (product[1] & (lshift(one, i-128)))) {
					product[2] ^= rshift(pp, 64-(i-128));
					product[3] ^= lshift(pp, i-128);
				}
			}
			ur2[j] = product[2];
			ur2[j+1] = product[3];
		}
	} else {
		for (j = 0; j < nbytes; j+=2) {
			product[0] = 0;
			product[1] = 0;
			product[2] = 0;
			product[3] = 0;
			for (i = 0; i < num_ones; i++) {
				if (scratch[i]<64) {
					product[1] ^= rshift(ur1[j], 64-scratch[i]);
					product[2] ^= lshift(ur1[j], scratch[i]);
					product[2] ^= rshift(ur1[j+1], 64-scratch[i]);
					product[3] ^= lshift(ur1[j+1], scratch[i]);
				} else if (scratch[i] >= 64) {
					product[0] ^= rshift(ur1[j], 128-scratch[i]);
					product[1] ^= lshift(ur1[j], scratch[i]-64);
					product[1] ^= rshift(ur1[j+1], 128-scratch[i]);
					product[2] ^= lshift(ur1[j+1], scratch[i]-64);
				}
			}
			for (i = top+128; i >= 128; i--) {
				if ((i >= 192) && (product[0] & (lshift(one, i-192)))) {
					product[1] ^= rshift(pp, 64-(i-192));
					product[2] ^= lshift(pp, i-192);
				} else if ((i >= 128) && (product[1] & (lshift(one, i-128)))) {
					product[2] ^= rshift(pp, 64-(i-128));
					product[3] ^= lshift(pp, i-128);
				}
			}
			ur2[j] ^= product[2];
			ur2[j+1] ^= product[3];
		}
	}
	return;

}

void galois_w128_bitmatrix_region_multiply (unsigned char *region,
		uint64_t multby[2],
		int nbytes,
		unsigned char *r2,
		int add)
{
	uint64_t *ur1, *ur2, *ur2top;
	uint64_t pp, one, prod[2], wshift, savem[2];
	int i, j;
	uint64_t scratch[128][2];

	ur1 = (uint64_t *) region;
	ur2 = (r2 == NULL) ? ur1 : (uint64_t *) r2;
	nbytes /= 8;
	ur2top = ur2 + nbytes;

	if (multby == 0) {
		if (!add) {
			while (ur2 < ur2top) {
				*ur2 = 0;
				ur2++;
			}
		}
		return;
	}

	one = 1;
	wshift = lshift(one, 63);
	pp = 0x87;
	
	savem[0] = multby[0];
	savem[1] = multby[1];

	for (i = 0; i < 128; i++) {
		scratch[i][0] = savem[0];
		scratch[i][1] = savem[1];

		if (savem[0] & wshift) {
			savem[0] = savem[0] << 1;
			if ((savem[1] & wshift) != 0) savem[0] = savem[0] ^ 1;
			savem[1] = savem[1] << 1;
			savem[1] = savem[1] ^ pp;
		} else {
			savem[0] = savem[0] << 1;
			if ((savem[1] & wshift) != 0) savem[0] = savem[0] ^ 1;
			savem[1] = savem[1] << 1;
		}
	}

	if (!add) {
		for (j = 0; j < nbytes; j+=2) {
			prod[0] = 0;
			prod[1] = 0;
			for (i = 0; i < 64; i++) {
				if ((ur1[j+1] & lshift(one, i))) {
					prod[0] ^= scratch[i][0];
					prod[1] ^= scratch[i][1];
				}
			}
			for (i = 0; i < 64; i++) {
				if ((ur1[j] & lshift(one, i))) {
					prod[0] ^= scratch[i+64][0];
					prod[1] ^= scratch[i+64][1];
				}
			}
			ur2[j] = prod[0];
			ur2[j+1] = prod[1];
		}
	} else {
		for (j = 0; j < nbytes; j+=2) {
			prod[0] = 0;
			prod[1] = 0;
			for (i = 0; i < 64; i++) {
				if ((ur1[j+1] & lshift(one, i))) {
					prod[0] ^= scratch[i][0];
					prod[1] ^= scratch[i][1];
				}
				if ((ur1[j] & lshift(one, i))) {
					prod[0] ^= scratch[i+64][0];
					prod[1] ^= scratch[i+64][1];
				}
			}
			ur2[j] ^= prod[0];
			ur2[j+1] ^= prod[1];
		}	
	}
	return;

}

void galois_w128_split8_table_region_multiply (unsigned char *region,
		uint64_t multby[2],
		int nbytes,
		unsigned char *r2,
		int add)
{
	uint64_t *ur1, *ur2;
	uint64_t *lp2, *lptop;
	uint64_t i, j, a, b, accumulator[2], i8, j8, k;
	uint64_t acache[16];

	ur1 = (uint64_t *) region;
	ur2 = (r2 == NULL) ? ur1 : (uint64_t *) r2;
	nbytes /= sizeof(uint64_t);

	if (multby == 0) {
		if (!add) {
			lp2 = (uint64_t *) ur2;
			ur2 += nbytes;
			lptop = (uint64_t *) ur2;
			while (lp2 < lptop) { *lp2 = 0; lp2++; }
		}
		return;
	}

	if (galois_w128_split8_tables[0]== NULL) {
		if (galois_w128_create_split8_tables() < 0) {
			fprintf(stderr, "galois_w128_split8_table_region_multiply -- couldn't make split multiplication tables\n");
			exit(1);
		}
	}

	/* If we're overwriting r2, then we can't do better than just calling split_multiply.
		 We'll inline it here to save on the procedure call overhead */

	i8 = 0;
	for (i = 0; i < 8; i++) {
		acache[i] = ((multby[1] >> i8) & 255) << 9;
		acache[i+8] = ((multby[0] >> i8) & 255) << 9;
		i8 += 8;
	}
	if (!add) {
		for (k = 0; k < nbytes; k+=2) {
			accumulator[0] = 0;
			accumulator[1] = 0;
			for (i = 0; i < 16; i++) {
				a = acache[i];
				j8 = 0;
				for (j = 0; j < 16; j++) {
					if (j < 8) b = ((ur1[k+1] >> j8) & 255) << 1;
					else b = ((ur1[k] >> (j8-64)) & 255) << 1;
					accumulator[0] ^= galois_w128_split8_tables[i+j][a|b];
					accumulator[1] ^= galois_w128_split8_tables[i+j][(a|b)+1];
					j8 += 8;
				}
			}
			ur2[k] = accumulator[0];
			ur2[k+1] = accumulator[1];
		}
	} else {
		for (k = 0; k < nbytes; k+=2) {
			accumulator[0] = 0;
			accumulator[1] = 0;
			for (i = 0; i < 16; i++) {
				a = acache[i];
				j8 = 0;
				for (j = 0; j < 16; j++) {
					if (j < 8) b = ((ur1[k+1] >> j8) & 255) << 1;
					else b = ((ur1[k] >> (j8-64)) & 255) << 1;
					accumulator[0] ^= galois_w128_split8_tables[i+j][a|b];
					accumulator[1] ^= galois_w128_split8_tables[i+j][(a|b)+1];
					j8 += 8;
				}
			}
			ur2[k] = accumulator[0];
			ur2[k+1] = accumulator[1];
		}
	}
	return;

}

uint64_t **galois_w128_create_lazy_split8_row(uint64_t multby[2]) {
	uint64_t p2, index, ishift, jshift, i, j, i8, j8;
	uint64_t p1elt[2], p2elt[2], acache[16], **table_row, *tmp;

	table_row = (uint64_t **) malloc(sizeof(uint64_t *) * 256);
	if (table_row == NULL) { return NULL; }

	for (i = 0; i < 256; i++) {
		table_row[i] = (uint64_t *) malloc(sizeof(uint64_t) * (1 << 8) * 2);
		if (table_row[i] == NULL) {
			for (i--; i >= 0; i--) free(table_row[i]);
			free(table_row);
			return NULL;
		}
	}

	i8 = 0;
	for (i = 0; i < 8; i++) {
		acache[i] = (multby[1] >> i8) & 255;
		acache[i+8] = (multby[0] >> i8) & 255;
		i8 += 8;
	}

	for (i = 0; i < 16; i++) {
		for (j = 0; j < 16; j++) {
			if (i+j > 15) {
				jshift = (i+j-15) * 8;
				ishift = 15 * 8;
			} else {
				jshift = (i+j) * 8;
				ishift = 0;
			}

			//p1elt = (acache[i] << ishift);
			index = 0;
			if (ishift != 0) {
				p1elt[0] = lshift(acache[i], ishift-64);
				p1elt[1] = 0;
			} else {
				p1elt[0] = 0;
				p1elt[1] = acache[i];
			}
			
			for (p2 = 0; p2 < 256; p2++) {
				if (jshift < 64) {
					p2elt[0] = 0;
					p2elt[1] = lshift(p2, jshift);
				} else {
					p2elt[0] = lshift(p2, jshift-64);
					p2elt[1] = 0;
				}
				tmp = galois_w128_bitmatrix_multiply(p1elt, p2elt);
				table_row[i*16+j][index] = tmp[0];
				table_row[i*16+j][index+1] = tmp[1];
				free(tmp);
				index+=2;
			}
		}
	}
	return table_row;
}

void galois_w128_lazy_split8_table_region_multiply (unsigned char *region,
		uint64_t multby[2],
		int nbytes,
		unsigned char *r2,
		int add)
{
	uint64_t *ur1, *ur2, *cp, *ur2top;
	uint64_t *lp2, *lptop;
	uint64_t i, j, b, accumulator[2], i8, j8, k;
	uint64_t **table_row;

	ur1 = (uint64_t *) region;
	ur2 = (r2 == NULL) ? ur1 : (uint64_t *) r2;
	nbytes /= sizeof(uint64_t);
	ur2top = ur2 + nbytes;

	if (multby == 0) {
		if (!add) {
			bzero(r2, nbytes * sizeof(uint64_t));
		}
		return;
	}

	/* If we're overwriting r2, then we can't do better than just calling split_multiply.
		 We'll inline it here to save on the procedure call overhead */

	table_row = galois_w128_create_lazy_split8_row (multby);
	i8 = 0;
	if (!add) {
		for (k = 0; k < nbytes; k+=2) {
			accumulator[0] = 0;
			accumulator[1] = 0;
			for (i = 0; i < 16; i++) {
				j8 = 0;
				for (j = 0; j < 16; j++) {
					if (j < 8) b = ((rshift(ur1[k+1], j8)) & 255) << 1;
					else b = ((rshift(ur1[k], j8-64)) & 255) << 1;
					
					accumulator[0] ^= table_row[i*16 + j][b];
					accumulator[1] ^= table_row[i*16 + j][b+1];
					j8 += 8;
				}
			}
			ur2[k] = accumulator[0];
			ur2[k+1] = accumulator[1];
		}
	} else {
		for (k = 0; k < nbytes; k+=2) {
			accumulator[0] = 0;
			accumulator[1] = 0;
			for (i = 0; i < 16; i++) {
				j8 = 0;
				for (j = 0; j < 16; j++) {
					if (j < 8) b = ((rshift(ur1[k+1], j8)) & 255) << 1;
					else b = ((rshift(ur1[k], j8-64)) & 255) << 1;
					
					accumulator[0] ^= table_row[i*16 + j][b];
					accumulator[1] ^= table_row[i*16 + j][b+1];
					j8 += 8;
				}
			}
			ur2[k]   ^= accumulator[0];
			ur2[k+1] ^= accumulator[1];
		}
	}
	for (i = 0; i < 256; i++) free(table_row[i]);
	free(table_row);
	return;
	
}

uint64_t **galois_w128_create_lazy_split16_row(uint64_t multby[2]) {
	int index, ishift, jshift, i, j, i16;
	uint64_t p2, p1elt[2], p2elt[2], acache[8], **table_row, *tmp;

	table_row = (uint64_t **) malloc(sizeof(uint64_t *) * 64);
	if (table_row == NULL) { return NULL; }

	for (i = 0; i < 64; i++) {
		table_row[i] = (uint64_t *) malloc(sizeof(uint64_t) * (1 << 16) * 2);
		if (table_row[i] == NULL) {
			for (i--; i >= 0; i--) free(table_row[i]);
			free(table_row);
			return NULL;
		}
	}

	i16 = 0;
	for (i = 0; i < 4; i++) {
		acache[i] = (rshift(multby[1], i16) & 65535);
		acache[i+4] = (rshift(multby[0], i16) & 65535);
		i16 += 16;
	}

	for (i = 0; i < 8; i++) {
		for (j = 0; j < 8; j++) {
			if (i+j > 7) {
				jshift = (i+j-7) * 16;
				ishift = 7 * 16;
			} else {
				jshift = (i+j) * 16;
				ishift = 0;
			}

			index = 0;
			if (ishift != 0) {
				p1elt[0] = lshift(acache[i], ishift-64);
				p1elt[1] = 0;
			} else {
				p1elt[0] = 0;
				p1elt[1] = acache[i];
			}

			for (p2 = 0; p2 < 65536; p2++) {
				if (jshift < 64) {
					p2elt[0] = 0;
					p2elt[1] = lshift(p2, jshift);
				} else {
					p2elt[0] = lshift(p2, jshift-64);
					p2elt[1] = 0;
				}
				tmp = galois_w128_bitmatrix_multiply(p1elt, p2elt);
				table_row[i*8+j][index] = tmp[0];
				table_row[i*8+j][index+1] = tmp[1];
				index+=2;
				free(tmp);
			}
		}
	}
	return table_row;
}

void galois_w128_lazy_split16_table_region_multiply (unsigned char *region,
		uint64_t multby[2],
		int nbytes,
		unsigned char *r2,
		int add)
{
	uint64_t *ur1, *ur2, *cp, *ur2top;
	uint64_t *lp2, *lptop;
	uint64_t i, j, b, accumulator[2], j16, k;
	uint64_t **table_row;

	ur1 = (uint64_t *) region;
	ur2 = (r2 == NULL) ? ur1 : (uint64_t *) r2;
	nbytes /= sizeof(uint64_t);
	ur2top = ur2 + nbytes;

	if (multby == 0) {
		if (!add) {
			bzero(r2, nbytes * sizeof(uint64_t));
		}
		return;
	}

	/* If we're overwriting r2, then we can't do better than just calling split_multiply.
		 We'll inline it here to save on the procedure call overhead */

	table_row = galois_w128_create_lazy_split16_row (multby);
	if (table_row == NULL) {
		fprintf(stderr, "galois_w128_lazy_split16_table_region_multiply -- couldn't make multiplication tables");
		exit(1);
	}
	if (!add) {
		for (k = 0; k < nbytes; k+=2) {
			accumulator[0] = 0;
			accumulator[1] = 0;
			for (i = 0; i < 8; i++) {
				j16 = 0;
				for (j = 0; j < 8; j++) {
					if (j < 4) b = (rshift(ur1[k+1], j16) & 65535) << 1;
					else b = (rshift(ur1[k], j16-64) & 65535) << 1;
					
					accumulator[0] ^= table_row[i*8 + j][b];
					accumulator[1] ^= table_row[i*8 + j][b+1];
					j16 += 16;
				}
			}
			ur2[k] = accumulator[0];
			ur2[k+1] = accumulator[1];
		}
	} else {
		for (k = 0; k < nbytes; k+=2) {
			accumulator[0] = 0;
			accumulator[1] = 0;
			for (i = 0; i < 8; i++) {
				j16 = 0;
				for (j = 0; j < 8; j++) {
					if (j < 4) b = (rshift(ur1[k+1], j16) & 65535) << 1;
					else b = (rshift(ur1[k], j16-64) & 65535) << 1;
					
					accumulator[0] ^= table_row[i*8 + j][b];
					accumulator[1] ^= table_row[i*8 + j][b+1];
					j16 += 16;
				}
			}
			ur2[k] ^= accumulator[0];
			ur2[k+1] ^= accumulator[1];
		}
	}
	for (i = 0; i < 64; i++) free(table_row[i]);
	free(table_row);
	return;

}

void galois_w128_composite_k2_region_multiply (unsigned char *region,
		uint64_t multby[2],
		int nbytes,
		unsigned char *r2,
		int add)
{
	uint64_t *ur1, *ur2, *lp, *lp2;
	int i, l;

	uint64_t prod[2];
	uint64_t s, p2, p3;
	uint64_t one;

	ur1 = (uint64_t *) region;
	ur2 = (r2 == NULL) ? ur1 : (uint64_t *) r2;
	nbytes /= 8;

	if (multby == 0) {
		if (!add) {
			lp = (uint64_t *) ur2;
			ur2 += nbytes;
			lp2 = (uint64_t *) ur2;
			while (lp < lp2) { *lp = 0; lp++; }
		}
		return;
	}

	if (galois_w64_split8_tables[0] == NULL) {
		if (galois_w64_create_split8_tables() < 0) {
			fprintf(stderr, "galois_w128_composite_k2_region_multiply: Error creating tables for w=64\n");
			exit(1);
		}
	}

	s = 2;
	l = 64;
	one = 1;

	if (r2 == NULL || !add) {
		for (i = 0; i < nbytes; i+=2) {

			p2 = galois_w64_split8_table_multiply(multby[0], ur1[i]);
			p3 = p2;
			p2 = galois_w64_split8_table_multiply(s, p2);

			prod[0]  = galois_w64_split8_table_multiply(multby[0], ur1[i+1]);
			prod[0] ^= galois_w64_split8_table_multiply(multby[1], ur1[i]);
			prod[0] ^= p2;
			prod[1]  = p3;
			prod[1] ^= galois_w64_split8_table_multiply(multby[1], ur1[i+1]);
			
			ur2[i]   = prod[0];
			ur2[i+1] = prod[1];
		}
	} else {
		for (i = 0; i < nbytes; i+=2) {

			p2 = galois_w64_split8_table_multiply(multby[0], ur1[i]);
			p3 = p2;
			p2 = galois_w64_split8_table_multiply(s, p2);

			prod[0]  = galois_w64_split8_table_multiply(multby[0], ur1[i+1]);
			prod[0] ^= galois_w64_split8_table_multiply(multby[1], ur1[i]);
			prod[0] ^= p2;
			prod[1]  = p3;
			prod[1] ^= galois_w64_split8_table_multiply(multby[1], ur1[i+1]);
			
			ur2[i]   ^= prod[0];
			ur2[i+1] ^= prod[1];
		}
	}
	return;
	
}

void galois_w128_split_composite_k2_region_multiply (unsigned char *region,
		uint64_t multby[2],
		int nbytes,
		unsigned char *r2,
		int add)
{
	/* low region pointer, high region pointer for destination and source */
	/* multby split into high, low, and s * high */
	int i, j;
	unsigned char *dhp, *shp, *dlp, *slp;
	uint64_t m1, m0, sm1;
	uint64_t one;
	unsigned char *tmp;

	/* set region pointers */
	shp = region;
	slp = region + nbytes/2;
	dhp = (r2 == NULL) ? shp : r2;
	dlp = dhp + nbytes/2;
	nbytes = nbytes/2;

	/* is multby 0? Speed that up */
	if (multby == 0) {
		if (!add) {
			bzero (dhp, nbytes);
			bzero (dlp, nbytes);
		}
		return;
	}

	/* create w64 split tables */
	if (galois_w64_split8_tables[0] == NULL) {
		if (galois_w64_create_split8_tables() < 0) {
			fprintf(stderr, "galois_w128_split_composite_k2_region_multiply -- couldn't make split tables\n");
			exit(1);
		}
	}

	/* set up constant multipliers */
	one = 1;
	m1 = multby[0];
	m0 = multby[1];

	sm1 = galois_w64_split8_table_multiply(m1, 2);

	if (r2 == NULL) {
		tmp = (unsigned char *) malloc(sizeof(unsigned char) * nbytes*2);
		dhp = tmp;
		dlp = tmp + nbytes;
		galois_w64_split8_table_region_multiply(slp, m1, nbytes, dhp, add);
		galois_w64_split8_table_region_multiply(shp, m0, nbytes, dhp, 1);
		galois_w64_split8_table_region_multiply(shp, sm1, nbytes, dhp, 1);

		galois_w64_split8_table_region_multiply(shp, m1, nbytes, dlp, add);
		galois_w64_split8_table_region_multiply(slp, m0, nbytes, dlp, 1);

		free (region);
		region = tmp;

		return;
	}

	/* w64 region multiplies */
	galois_w64_split8_table_region_multiply(slp, m1, nbytes, dhp, add);
	galois_w64_split8_table_region_multiply(shp, m0, nbytes, dhp, 1);
	galois_w64_split8_table_region_multiply(shp, sm1, nbytes, dhp, 1);

	galois_w64_split8_table_region_multiply(shp, m1, nbytes, dlp, add);
	galois_w64_split8_table_region_multiply(slp, m0, nbytes, dlp, 1);

	return;

}

void galois_w128_recursive_composite_k2_region_multiply (unsigned char *region,
		uint64_t multby[2],
		int nbytes,
		unsigned char *r2,
		int add)
{
	/* low region pointer, high region pointer for destination and source */
	/* multby split into high, low, and s * high */
	int i, j;
	unsigned char *dhp, *shp, *dlp, *slp;
	uint64_t m1, m0, sm1;
	unsigned char *tmp;

	/* set region pointers */
	shp = region;
	slp = region + nbytes/2;
	dhp = (r2 == NULL) ? shp : r2;
	dlp = dhp + nbytes/2;
	nbytes = nbytes/2;

	/* is multby 0? Speed that up */
	if (multby[0] == 0 & multby[1] == 0) {
		if (!add) {
			bzero (dhp, nbytes);
			bzero (dlp, nbytes);
		}
		return;
	}

	/* set up constant multipliers */
	m1 = multby[0];
	m0 = multby[1];

	sm1 = galois_w64_recursive_composite_k2_multiply(m1, 2);

	if (r2 == NULL) {
		tmp = (unsigned char *) malloc(sizeof(unsigned char) * nbytes*2);
		dhp = tmp;
		dlp = tmp + nbytes;
		galois_w64_recursive_composite_k2_region_multiply(slp, m1, nbytes, dhp, add);
		galois_w64_recursive_composite_k2_region_multiply(shp, m0, nbytes, dhp, 1);
		galois_w64_recursive_composite_k2_region_multiply(shp, sm1, nbytes, dhp, 1);

		galois_w64_recursive_composite_k2_region_multiply(shp, m1, nbytes, dlp, add);
		galois_w64_recursive_composite_k2_region_multiply(slp, m0, nbytes, dlp, 1);

		free (region);
		region = tmp;

		return;
	}

	/* w64 region multiplies */
	galois_w64_recursive_composite_k2_region_multiply(slp, m1, nbytes, dhp, add);
	galois_w64_recursive_composite_k2_region_multiply(shp, m0, nbytes, dhp, 1);
	galois_w64_recursive_composite_k2_region_multiply(shp, sm1, nbytes, dhp, 1);

	galois_w64_recursive_composite_k2_region_multiply(shp, m1, nbytes, dlp, add);
	galois_w64_recursive_composite_k2_region_multiply(slp, m0, nbytes, dlp, 1);

	return;

}

/* Functions that do not directly depend upon another galois implementation */

void galois_region_xor(unsigned char *r1,         /* Region 1 */
		unsigned char *r2,         				  /* Region 2 */
		unsigned char *r3,         				  /* Sum region (r3 = r1 ^ r2) -- can be r1 or r2 */
		int nbytes)               				  /* Number of bytes in region */ 
{
	uint64_t *l1;
	uint64_t *l2;
	uint64_t *l3;
	uint64_t *ltop;
	unsigned char *ctop;

	ctop = r1 + nbytes;
	ltop = (uint64_t *) ctop;
	l1 = (uint64_t *) r1;
	l2 = (uint64_t *) r2;
	l3 = (uint64_t *) r3;

	while (l1 < ltop) {
		*l3 = ((*l1)  ^ (*l2));
		l1++;
		l2++;
		l3++;
	}
}

void reed_sol_galois_w08_region_multby_2(unsigned char *region, int nbytes) {
  uint32_t *l1;
  uint32_t *ltop;
  unsigned char *ctop;
  uint32_t tmp, tmp2;
  int prim08 = -1;
  int mask08_1 = -1;
  int mask08_2 = -1;

  if (prim08 == -1) {
    tmp = galois_single_multiply((1<<7), 2, 8);
    prim08 = 0;
    while (tmp != 0) {
      prim08 |= tmp;
      tmp = (tmp << 8);
    }
    tmp = (1 << 8) - 2;
    mask08_1 = 0;
    while (tmp != 0) {
      mask08_1 |= tmp;
      tmp = (tmp << 8);
    }
    tmp = (1 << 7);
    mask08_2 = 0;
    while (tmp != 0) {
      mask08_2 |= tmp;
      tmp = (tmp << 8);
    }
  }

  ctop = region + nbytes;
  ltop = (uint32_t *) ctop;
  l1 = (uint32_t *) region;

  while (l1 < ltop) {
    tmp = ((*l1) << 1) & mask08_1;
    tmp2 = (*l1) & mask08_2;
    tmp2 = ((tmp2 << 1) - (tmp2 >> 7));
    *l1 = (tmp ^ (tmp2 & prim08));
    l1++;
  }
}

void reed_sol_galois_w16_region_multby_2(unsigned char *region, int nbytes) {
  uint32_t *l1;
  uint32_t *ltop;
  unsigned char *ctop;
  uint32_t tmp, tmp2;

  int prim16 = -1;
  int mask16_1 = -1;
  int mask16_2 = -1;


  if (prim16 == -1) {
    tmp = galois_single_multiply((1 << 15), 2, 16);

    prim16 = 0;
    while (tmp != 0) {
      prim16 |= tmp;
      tmp = (tmp << 16);
    }
    tmp = (1 << 16) - 2;
    mask16_1 = 0;
    while (tmp != 0) {
      mask16_1 |= tmp;
      tmp = (tmp << 16);
    }
    tmp = (1 << 15);
    mask16_2 = 0;
    while (tmp != 0) {
      mask16_2 |= tmp;
      tmp = (tmp << 16);
    }
  }

  ctop = region + nbytes;
  ltop = (uint32_t *) ctop;
  l1 = (uint32_t *) region;

  while (l1 < ltop) {
    tmp = ((*l1) << 1) & mask16_1;
    tmp2 = (*l1) & mask16_2;
    tmp2 = ((tmp2 << 1) - (tmp2 >> 15));
    *l1 = (tmp ^ (tmp2 & prim16));
    l1++;
  }
}

void reed_sol_galois_w32_region_multby_2(unsigned char *region, int nbytes) {
  uint32_t *l1;
  uint32_t *ltop;
  unsigned char *ctop;
  int prim32 = -1;

  if (prim32 == -1) prim32 = galois_single_multiply((1 << 31), 2, 32);

  ctop = region + nbytes;
  ltop = (uint32_t *) ctop;
  l1 = (uint32_t *) region;

  while (l1 < ltop) {
    *l1 = ((*l1) << 1) ^ ((*l1 & 0x80000000) ? prim32 : 0);
    l1++;
  }
}
